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ABSTRACT 


This  thesis  describes  an  analytical  study  of  in-situ  creep  behavior  in 
ice-rich  permafrost.  An  incremental  initial  strain  procedure  was  formulated 
to  solve  steady  state  creep  problems  using  the  finite  element  method.  Two 
case  histories  were  analyzed  where  naturally  occurring  creep  had  been 
monitored  for  several  years. 

The  numerical  analyses  of  the  left  bank  of  the  Great  Bear  River  at  the 
proposed  Arctic  Gas  crossing  has  shown  that  the  steady  state  creep  in  the 
ice-rich  g I ac i o I acustr i ne  clay  at  the  site  can  be  modelled  by  a  simple  power 
law  with  an  exponent  of  3.0  and  a  coefficient  of  3.33  x  1 0-9  kPa-^  yr_1  . 
This  strain  rate  is  six  times  slower  than  the  value  for  polycrystalline  ice 
at  an  equivalent  temperature.  The  exact  form  of  the  constitutive 
relationship  for  the  g I ac iode I ta i c  sand  overlying  the  clay  remains  unclear. 

A  review  of  the  available  in-situ  deformation  studies  carried  out  at 
the  Fox  Tunnel  near  Fairbanks,  Alaska,  showed  that  attenuating  creep 
behavior  extended  well  over  one  year  and  contributed  significantly  to  the 
total  room  closure.  In  this  case,  the  flow  law  for  polycrystalline  ice  did 
not  yield  an  upper  bound  solution  to  the  observed  room  closure  measurements. 
However,  it  is  felt  that  the  difference  lies  in  time-dependent  failure 
caused  by  stress  release  which  was  not  accounted  for  in  the  numerical 
mode  I  . 
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CHAPTER  I 


INTRODUCTION 


1 . 1  Genera  I 


Depletion  of  conventional  oil  and  gas  reserves  has  led  to  an  increase 
in  resource  exploration  activity  in  the  northern  frontier  lands  of  Canada 
and  the  United  States  during  the  past  decade  and  into  the  foreseeable 
future.  Rising  world  energy  prices  has  enhanced  the  economic  feasibility  of 
exploiting  these  reserves  and  transporting  them  to  market.  The  largest  such 
project  completed  to  date  is  the  Alyeska  Pipeline  which  transports  Alaskan 
North  Slope  oil  to  a  shipping  terminal  at  Valdez.  The  first  major  northern 
oil  pipeline  to  be  completed  in  Canada  will  be  the  900  km  long  pipeline 
forecast  to  be  in  place  by  1985  to  transport  oil  south  from  the  expanded  oil 
field  in  Norman  Wells,  N . W.T.  to  Zama,  Alberta  where  it  will  connect  with  an 
existing  p ipe I i ne. 

Resource  exploration  and  development  will  require  major  civil  works  to 
be  constructed  on  or  within  permafrost.  These  structures  will  present 
unique  problems  to  engineers  as  the  strength  and  deformation  behavior  in 
frozen  soil  will  respond  in  both  the  time  and  thermal  domain.  The  presence 
of  ice  as  either  discrete  segregated  structures  or  in  the  pore  spaces, 
imparts  a  major  influence  on  the  creep  behavior  that  character i zes  frozen 
soil. 


The  i  n-s  i  tu  creep  behavior  of  frozen  soil  will  have  a  major  influence 
on  establishing  geotechnical  design  guidelines  that  ensure  serviceable 
foundation  performance  for  structures  founded  in  permafrost.  Ladanyi  (1972) 
has  developed  an  engineering  theory  for  creep  in  frozen  soils.  The  finite 
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element  technique  is  now  routinely  used  for  solving  a  large  class  of 
engineering  problems.  The  aim  of  the  research  reported  herein  was  to  use 
the  finite  element  technique  to  predict  the  in-situ  creep  deformation 
patterns  in  ice-rich  permafrost  foundation  soils. 

To  date,  these  are  only  two  documented  case  histories  comparing  long 
term  in-situ  and  I aboratory  creep  behavior  of  frozen  soil.  An  analytical 
study  was  carried  out  by  Thompson  and  Sayles  (1972)  in  conjunction  with  a 
field  program  of  the  U.S.  Army  Cold  Regions  Research  and  Engineering 
Laboratory  (USA  CRREL)  in  the  Fox  Tunnel  north  of  Fairbanks,  Alaska.  In- 
situ  measurement  of  creep  deformations  were  carried  out  for  three  and  one 
half  years.  Savigny  (1980)  reported  the  results  of  an  extensive  field  and 
laboratory  program  on  the  proposed  Arctic  Gas  crossing  of  the  Great  Bear 
River  near  Fort  Norman,  N.W.T.  Naturally  occurring  creep  movements  in  the 
left  bank  of  the  Great  Bear  River  were  recorded  for  a  two  year  period. 

1 .2  Scope  of  Thesis 

Chapter  I  I  presents  a  review  of  the  rheological  behavior  of  frozen 
soil.  Since  the  deformation  behavior  of  frozen  soil  is  governed  by  the 
amount  of  ice  present,  a  detailed  review  of  the  deformation  behavior  of 
po  I  ycrysta  II  i  ne  ice  is  presented.  This  is  followed  by  a  review  of  the 
deformation  behavior  for  reconstituted  frozen  soil  and  natural  permafrost. 
The  existing  constitutive  relationships  for  creep  in  natural  permafrost 
samples  are  examined.  An  upper  bound  flow  law  for  the  creep  in  ice-rich 
fine-grained  permafrost  is  presented  at  the  conclusion  of  this  chapter. 

Chapter  III  presents  the  theoretical  basis  for  the  incremental  initial 
strain  procedure  for  solving  for  the  stress  state  in  a  non-linear  viscous 
medium.  The  treatment  is  restricted  to  steady  state  creep  problems  since 
the  empirically  developed  constitutive  relationship  for  frozen  soil  is  given 
by  a  simple  power  law  relating  strain  rate  to  stress.  The  development  of 
the  plane  strain  finite  element  programme  CREEP  is  also  presented. 
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Finite  element  analyses  of  the  in-situ  creep  behavior  of  the  left  bank 
of  the  Great  Bear  River  at  the  proposed  Arctic  Gas  crossing  is  presented  in 
Chapter  IV.  These  analyses  were  carried  out  in  an  attempt  to  verify  the 
simple  power  law  as  the  constitutive  equation  governing  the  behavior  of 
frozen  ice-rich,  fine-grained  permafrost  soils. 

The  in-situ  deformation  behavior  at  the  Fox  Tunnel  at  Fairbanks,  Alaska 
is  examined  in  Chapter  V.  A  synthesis  of  al I  the  in-situ  deformation 
studies  reported  in  the  literature  for  the  Fox  Tunnel  was  carried  out  in  an 
attempt  to  gain  some  insight  into  the  deformation  processes  surrounding 
underground  openings  in  frozen  ground. 

The  final  chapter  presents  a  summary  of  the  results  of  this  research  as 
well  as  the  limitations  of  the  current  analytical  approach.  Recommended 
suggestions  for  future  analytical  and  laboratory  creep  studies  in  frozen 
soils  are  also  presented. 

The  formulation  of  the  finite  element  equations  for  solving  creep 
problems  with  the  incremental  initial  strain  procedure  is  presented  in 
Appendices  A  and  B.  A  user’s  manual  and  programme  listing  for  the  finite 
element  programme  CREEP  are  given  in  Appendices  C  and  D,  respectively.  The 
accuracy  of  the  finite  element  programme  is  compared  with  a  thick  wall 
cylinder  closed  form  solution  in  Appendix  E. 
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CHAPTER  I  I 


RHEOLOGICAL  BEHAVIOR  OF  FROZEN  SOILS 


2 . 1  Genera  I 

An  understanding  of  the  deformation  behavior  of  frozen  ground  over  a 
wide  range  of  stress  and  temperature  conditions  is  reguired  before  a 
numerical  analysis  can  be  carried  out.  Deformation  behavior  of  frozen 
ground  is  generally  governed  by  the  amount  of  ice  present  and  its 
temperature . 

Laboratory  studies  over  the  past  two  decades  have  focussed  attention  on 
deriving  empirical  constitutive  relationships  describing  the  stress-stra i n- 
time  behavior  of  frozen  soils.  Many  of  the  early  studies  were  restricted  to 
short  duration  creep  tests  on  ice  and  remoulded  frozen  soils  at  stresses  and 
temperatures  well  beyond  the  range  of  interest  of  practical  geotechnical 
problems.  As  the  major  factors  influencing  the  creep  behavior  of  frozen 
soil  have  been  delineated,  specifically  the  soil  temperature  and  ice 
structure,  the  need  for  high  guality,  long  term,  low  stress  creep  tests  on 
natural  permafrost  samples  has  arisen. 

The  following  sections  will  present  a  brief  review  of  the  time- 
dependent  load-deformation  behavior  of  ice,  remoulded  frozen  soil  and 
natural  permafrost  soils  as  it  relates  to  low  stress  creep  behavior.  Unless 
otherwise  mentioned,  stress  and  temperature  conditions  commonly  encountered 
in  geotechnical  engineering  are  assumed. 
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2.2 


Composition  of  Frozen  Soil 


Frozen  soil  can  be  considered  as  a  complex,  four  phase,  natural 
formation,  consisting  of  solid  mineral  particles,  ice,  unfrozen  water  and 
gases.  The  solid  mineral  particles  present  in  a  frozen  soil  have  an 
important  influence  on  its  geotechnical  properties.  The  primary  factors  are 
the  grain  size  and  shape  as  well  as  the  physicochemical  nature  of  the 
mineral  surfaces,  which  are  determined  by  the  m i nera I og i ca I  composition  and 
the  cations  present.  The  gaseous  phase  consists  of  water  vapour  and  has  an 
insignificant  influence  on  the  behavior  of  frozen  soils. 


Unfrozen  water  is  present  in  two  states  as  strongly  bound  and  loosely 
bound  water.  The  strongly  bound  water  is  adjacent  to  the  soil  particle. 
The  very  high  e  I  ectro-mo  I  ecu  I  ar  forces  present  in  this  layer  of  water 

suppresses  the  formation  of  ice  crystals,  even  at  very  low  temperatures. 
The  i n termo I  ecu  I ar  forces  are  also  present  in  the  loosely  bound  water 

surrounding  the  strongly  bound  water.  However,  this  layer  is  capable  of 
releasing  the  heat  of  crystallization  at  temperatures  below  0°C. 

The  amount  of  unfrozen  water  present  in  a  frozen  soil  depends  on  its 
temperature,  specific  surface  and  type  of  soil  mineral  present  and  chemistry 
of  the  pore  water.  Each  frozen  soil  is  character  i  zed  by  a  specific  curve 
relating  its  unfrozen  water  content  and  temperature. 

Of  all  the  phases  present  in  frozen  soil,  ice  is  the  most  important 
component  in  the  rate,  time  and  temperature  dependent  properties  which 
characterize  frozen  soil.  Ice  can  be  present  both  as  discrete  segregated 
structures  and  in  the  pore  spaces  of  the  soil  mass.  Ice  present  in  frozen 
soils  is  generally  polycrystalline. 
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The  reader  is  referred  to  Anderson  and  Morgenstern  (1973)  for  a 
complete  summary  of  conditions  affecting  the  state  of  frozen  ground. 

2.3  Creep  of  Ice 


The  presence  of  ice  dominates  the  time-dependent  character i st i cs  of 
frozen  soil.  Ice  1H  is  the  predominant  ice  type  found  in  frozen  soil.  The 
structure  of  ice  1H  is  hexagonal.  Each  crystal  consists  of  layers  of 
hexagonal  rings,  referred  to  as  the  basal  planes.  The  perpend icu I ar 
direction  to  the  basal  plane  is  the  optic  or  1 C T  axis. 

Single  crystals  of  ice  deform  readily  at  low  stresses  along  discrete 
bands  parallel  to  the  basal  plane  of  the  crystal  structure.  This  type  of 
deformation  is  referred  to  as  basal  glide  (easy  glide).  Single  ice  crystals 
will  also  deform  in  the  perpend i cu I ar  or  non-basal  direction  (hard  glide). 
However,  for  a  given  strain  rate,  the  appl  ied  stress  can  be  as  much  as  20 
times  greater  for  hard  glide  than  for  easy  glide. 

Po I ycrysta I  I i ne  ice  is  composed  of  randomly  oriented  single  ice 
crystals.  Physical  processes  observed  during  the  deformation  of 
polycrystalline  ice  are  dislocation  climb,  grain  boundary  slip,  cavity 
formation  at  the  grain  boundaries,  recrysta I  I i zat ion  and  microcracks 
developing  into  ice  grains  (Michel,  1978). 

The  deformation  response  of  polycrystal  I ine  ice  subjected  to  a  load 
consists  of  four  distinct  regions  as  illustrated  by  curve  II  in  Figure  2.1. 
Each  of  these  regions  consist  of: 

1)  instantaneous  elasfic  strain  (0A) 

2)  creep  strain  in  a  primary  inode  at  a  decreasing  strain  (AB) 
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3)  secondary  creep  at  a  constant  minimum  strain  rate  (BC) 

4)  tertiary  creep  at  an  accelerating  strain  rate  leading  to 
fa i I ure  (CD) 

The  shape  of  the  creep  curve  is  a  function  of  the  stress  level  and  ice 
temperature.  For  a  given  temperature,  secondary  creep  is  suppressed  at  high 
stress  levels  (undamped  creep)  while  primary  creep  (damped  creep)  dominates 
at  low  stress  levels. 

Creep  tests  on  polycrystalline  ice  in  compression,  tension,  shear  and 
other  specialized  loading  configurations  have  been  carried  out  to  evaluate 
the  form  of  the  constitutive  equation.  Glen  (1975)  has  recognized  two 
possible  ways  of  identifying  a  flow  law.  The  first  flow  law  would  relate 
stress  to  strain  rate  at  a  very  long  time  under  the  application  of  load 
allowing  recrysta I  I i zat ion  to  occur.  The  second  form  of  the  flow  law 
relates  the  secondary  strain  rate  to  stress  before  any  recrysta II i zat ion  of 

the  ice.  At  higher  stresses,  the  first  flow  law  would  produce  deformation 
rates  greater  than  the  minimum  (or  secondary)  rate  because  recrysta II i zat i on 
may  produce  a  grain  size  and  orientation  more  favorable  for  plastic  flow. 

The  constitutive  equation  of  polycrystalline  ice  is  most  commonly 
represented  empirically  by  the  simple  power  law.  Experimental  evidence  has 
shown  that  the  power  creep  law  represents  the  steady  state  creep  data  in  the 
low  to  intermediate  stress  range.  From  a  practical  viewpoint,  this 
mathematical  representat ion  of  creep  behavior  is  simpler  to  use  than  a 
physical  theory  (i.e.  rate  process  theory)  because  the  material  parameters 
are  kept  to  a  minimum. 
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The  simple  power  law  relating  strain  rate  to  stress  is  given  by: 


e  =  A a  (2.1 ) 

where  £  =  axial  strain  rate  (time)-^ 

a  =  axial  stress 

A  =  coefficient  (function  of  temperature  and  ice  type) 
n  =  creep  exponent 

In  order  to  extend  uniaxial  creep  data  to  rnulti-axial  stress  states,  it 
is  common  practice  to  adopt  the  concept  of  effective  stress,  a  f  and 
effective  creep  strain  rate,  £e.  Several  definitions  of  effective  stress 
and  strain  rate  have  been  put  forward  in  the  metal  and  ice  creep  I iterature. 
Some  of  these  are  listed  in  Table  2.1. 

The  definition  of  effective  stress  is  defined  by  the  square  root  of  the 
second  invariant  of  the  deviatoric  stress  tensor  multiplied  by  a  constant. 
Likewise,  the  effective  strain  rate  is  defined  in  a  similar  manner  except 
that  the  second  invariant  of  the  strain  rate  tensor  is  used.  The  difference 
between  the  alternate  definitions  lies  in  the  value  of  the  constant.  Since 
the  stress  is  raised  to  an  exponent  in  the  simple  power  law,  the  magnitude 
of  the  constant  will  influence  the  effective  strain  rate  in  a  non-linear 
manner  as  shown  in  Figure  2.3.  As  shown  in  the  figure,  the  effective  strain 
rate  can  change  by  a  complete  order  of  magnitude  depending  on  which 
definition  of  effective  stress  is  chosen.  The  steady  state  creep  law  for 
ice  at  -2°C  ( Morgens tern  et  al.,  1980)  is  plotted  in  Figure  2.3. 
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TABLE  2.1  DEFINITIONS  OF  EFFECTIVE  STRESS  AND 
EFFECTIVE  STRAIN  RATE 


Source 

Effective  Stress 

Effective  Strain 

Rate 

Odqv i st  ( 1966 ) 

Dorn  et  al  (1945) 

Ladany i  (  1972 ) 

Finn ie  &  Heller  (1959) 

-2_Vj2 
^  ' 

Nye  (  1957) 

VI 

Meier  (1959) 

Vfc 

Vialov  ( 1962 ) 

VI 

vs  VI 

where 


a 

Jo  =  second  invariant  of  deviator ic  stress  tensor 

,e 

Jo  =  second  invariant  of  strain  rate  tensor 

Z. 
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In  this  work,  the  effective  stress  and  strain  rate  as  defined  by 
Odqvist  (1956)  has  been  used,  i.e.: 


(2.2) 


(2.3) 


These  definitions  are  adopted  because  they  recover  Norton's  law,  i.e. 
equation  2.1,  for  the  special  case  of  uniaxial  stress.  The  flow  law  now 
becomes : 


(2.4) 


Since  the  effective  stress  is  a  function  of  the  deviatoric  stress  only, 
the  strain  rate  is  independent  of  the  hydrostatic  state  of  stress. 

An  extensive  review  of  the  secondary  creep  data  of  polycrystalline  ice 
was  carried  out  by  Morgenstern  et  al.  (1980).  Based  on  this  review,  the 
authors  presented  the  values  shown  in  Table  2.2  for  the  parameters  A  and  n 
in  the  power  law. 


A  more  detailed  review  of  the  rheological  character i st i cs  of  polycrystalline 
ice  has  been  prepared  by  Glen  (1975),  Roggensack  (1977)  and  Sego  (1980). 


TABLE  2.2  CREEP  PARAMETERS  FOR  POLYCRYSTALLINE  ICE 


1 1 


Temperature  (°C) 

A( kPa  3yr  ) ) 

n 

-1 

CO 

i 

o 

X 

• 

3.0 

-2 

2.0  x  10"3 

3.0 

-5 

CO 

1 

o 

X 

o 

• 

3.0 

-10 

5.6  x  10"9 

3.0 

2.4  Creep  of  Frozen  Soi I 

Vialov  (1963)  described  the  physical  processes  involved  in  the  time 
dependent  plastic  deformation  of  frozen  soils.  The  application  of  a  load 
causes  local  stress  concentrations  at  soil-ice  contacts  leading  to  pressure 
melting  of  the  ice.  A  pressure  gradient  resulting  from  differences  in 
surface  tensions  causes  the  water  to  migrate  to  regions  of  lower  stress 
where  it  then  refreezes.  The  pressure  melting  process  is  accompanied  by  a 
breakdown  of  the  ice  and  structural  bonds  of  the  soil  grains  and  plastic 
deformation  of  the  ice,  causing  a  rearrangement  of  the  mineral  particles. 
The  end  result  is  the  process  known  as  creep. 

Many  laboratory  studies  have  been  carried  out  in  order  to  determine  the 
empirical  creep  parameters  for  reconsituted  frozen  soils  and  natural 
permafrost  soils.  The  creep  tests  carried  out  on  remoulded  frozen  soils  are 
helpful  in  determining  the  various  factors  that  influence  the  creep  rate  of 
frozen  soils.  However,  many  of  these  tests  were  carried  out  beyond  the 
temperature  and  stress  range  of  practical  interest  in  geotechnical 
eng i neer i ng. 
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A  second  difficulty  encountered  when  applying  the  creep  parameters  of 
reconstituted  frozen  soils  to  natural  permafrost  soils  is  reproducab i I i ty  of 
the  ice  facies  in  the  remoulded  samples.  In  natural  permafrost  samples,  the 
i ce  '5  PreSGr>t  in  the  pore  spaces  as  well  as  in  the  form  of  reticulate  and 
segregated  ice.  Laboratory  evidence  of  Savigny  (1980)  shows  deformations 
localized  along  segregated  ice  in  natural  permafrost  soils.  Since  the 
reticulate  ice  structure  can  not  be  reproduced  in  reconstituted  frozen 
samples,  natural  permafrost  samples  must  be  tested. 

Vialov  (1959)  describes  the  stress-stra i n  response  of  a  frozen  soil  as 
damped  or  undamped  creep.  Damped  deformations  occur  when  the  applied  stress 
is  less  than  the  long  term  strength  of  the  frozen  soil.  Damped  creep 
behavior  is  character i zed  by  an  instantaneous  elastic  displacement  followed 
by  primary  creep  where  the  deformation  rate  continuously  decreases  with 
line.  Curve  III  in  Figure  2.1  is  typical  of  damped  creep  behavior. 
Undamped  deformations  occur  when  the  applied  stress  exceeds  the  long  term 
strength.  In  undamped  creep,  the  transient  stage  gives  way  to  a  steady 
state  period  where  the  strain  rate  reaches  a  minimum  value.  The 
deformations  continue  to  grow  until  they  reach  a  limiting  value  or  the  onset 
of  tertiary  creep  where  the  deformation  rate  will  accelerate  to  failure.  A 

typical  undamped  deformation  response  is  given  by  curves  I  and  II  in  Figure 

2.1. 


Damped  creep  deformations  are  prevalent  in  dense,  ice-poor  frozen 
soils.  Since  ice  has  a  zero  long  term  strength,  its  deformation  behavior  is 
character i zed  by  undamped  creep.  The  deformation  response  of  ice-rich 
frozen  soils  will  lie  somewhere  between  these  two  limiting  cases. 


•  ft' » 


V 

IIS  t  •  •  ’►*  b:,i;fr*c  ^^1 
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Weaver  (1979)  proposed  a  classification  system  for  frozen  soils  based 
on  frozen  bulk  density  that  relates  to  low  stress  creep  behavior.  This 
classification  system  is  shown  in  Table  2.3. 

A  generalized  theory  of  creep  of  frozen  soil  put  forward  by  Ladanyi 
(1972)  advocates  a  mathematical  rather  than  a  physical  treatment  of  undamped 
creep  behavior.  The  proposed  uniaxial  power  creep  law  describing  steady 
state  creep  is  defined  as: 


e  -  e  (a  / a  ) 
e  c  e  c 


(2.5) 


where:  ac  -  temperature  dependent  creep  coefficient 

ec  =  secondary  strain  rate  corresponding 
to  the  application  of  the  stress,  0 

n  =  creep  exponent 


This  form  of  the  power  law  is  similar  to  that  given  by  equation  2.4  except 
that  if  is  given  in  a  normalized  form.  Making  use  of  the  definitions  of 
effective  stress  (equation  2.2),  the  strain  rate  for  an  axially  symmetric 
state  of  stress  ( i .e.  a 2  =  a 3)  is  given  by: 

ee  =  ^c[ (a1 “a3)^acI R  (2.6) 


a1  “  a3  = 


where: 


principal  stress  difference 


■ 


TABLE  2.3  FROZEN  SOIL  CLASSIFICATION  SYSTEM 
(Proposed  by  Weaver,  1979) 


SOIL  TYPE 


DESCR  IPT ION 


Dirty  Ice 


appl  ies  to  ice  that  has  a  low  sol  ids 

concentrat i on 

y  =  0.9  -  1.0  Mg/ 

the  soil  particles  present  reduce 

the  average  grain  size  of  the  ice 

crystals  resulting  in  higher  creep 

rates  than  pure  ice 


Very  Dirty  Ice 


Ice-Poor  Frozen 
So 


Ice-Rich  Frozen 
So 


.3 


applies  to  ice  that  has  medium  to 
high  solids  concentration 
y  ~  0.9  -  1  .0  and  1.6  -  1.8  Mg/rrr 
very  little  grain  to  grain  contact 
between  soil  particles 
lower  secondary  creep  rates  than 
polycrystalline  ice  because  soil 
mpedes  dislocation  movement 

applies  to  saturated  frozen  soil 
whose  deformation  patterns  are 
character i zed  by  primary  creep 
y  -  1 • 7  -  1.8  and  1.9  -  2.0  Mg/m^ 

applies  to  soils  that  have  a 
continuous  network  of  segregated  ice 
the  overall  creep  response  is 
complex  and  is  very  sensitive  to  the 
reticulate  structure  of  the  segregated 
ice,  bulk  density  grain  size  and 
ground  temperature 
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The  strain  rate  derived  from  equation  2.6  is  independent  of  the 
hydrostatic  pressure.  The  hydrostatic  stress  can  be  expected  to  influence 
the  stress-stra i n  rate  behavior  up  to  failure  of  unconsolidated  frictional 
earth  materials.  Ladanyi  (1972)  extended  the  stress-stra i n  rate 
constitutive  law  to  account  for  the  mean  normal  stress  using  a  two  or  three 
parameter  failure  theory. 

Vialov  (1962)  proposed  that  the  strength  of  a  frozen  soil  can  be 
represented  by  a  series  of  Mohr  failure  envelopes,  each  one  corresponding  to 
a  specific  time  to  failure.  Assuming  validity  of  the  Mohr-Coulomb  failure 
criterion  in  the  pre-failure  state,  the  dependence  of  strain  rate  on  mean 
normal  pressure  for  an  axially  symmetric  stress  state  is  expressed  as: 


e 

e 


e 

c 


( f+2) (a  ^-03) 
c 


3  (  f  —  1 )  a 


m 


n 


(2.7) 


where:  f  =  (1  +  sin<j>)/(1  -  sin<|>) 

am  =  (a1 +a2+CT3)/^ 

ij)  =  angle  of  internal  friction 

For  a  frictionless  soil,  f  =  1,  equation  2.7  simplifies  to  equation 
2.6.  Equation  2.7  assumes  full  mobilization  of  internal  friction  over  the 
pre-failure  state.  This  assumption  leads  to  a  non-zero  strain  rate  at  a 
zero  stress  difference.  Thus,  application  of  equation  2.7  should  be 
restricted  to  strains  close  to  failure.  This  limitation  may  also  be 
overcome  by  assuming  a  time  dependent  angle  of  internal  friction. 

Ladanyi  (1972)  used  the  extended  von  Mises  failure  criterion  to  account 
for  effect  of  mean  normal  stress  on  strain  rate  for  frozen  soils  with  low 
internal  friction.  For  an  axially  symmetric  state  of  stress,  the  strain 


■ 
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rate  can  be  expressed  in  terms  of  mean  normal  stress  and  principal  stress 
d  i  f f erence  as : 


e 

p 


n 


e 

c 


(r+1 ) (01-03)  -  3 ( r —  1 )om 
2  a 

c 


(2.8) 


where  r  =  ratio  of  uniaxial  compressive 

strength  to  uniaxial  tensile  strength 

Equation  2.8  is  subject  to  the  same  limitations  as  equation  2.7  in  that 
a  zero  stress  difference  gives  a  non-zero  strain  rate  assuming  a  constant 
strength  ratio,  r.  The  limitation  of  equations  2.7  and  2.8  are  overcome  by 

simply  assuming  the  strain  rate  in  a  pre-failure  state  is  independent  of 
hydrostatic  stress  and  only  the  creep  strength  is  made  a  function  of  the 

mean  normal  stress.  This  is  consistent  with  the  experimental  findings  of 
Sayles  (1973)  at  higher  strain  rates  and  confining  pressures. 

From  a  practical  viewpoint,  the  engineering  theory  is  more  favourable 
than  a  theory  describing  creep  in  terms  of  more  fundamental  laws  of  physics 
because  the  required  number  of  experimental  mater, ial  parameters  are  kept  to 
a  mini  mum. 

Tests  on  ice-rich  remoulded  frozen  soils  have  been  carried  by  Vialov 
(1962),  Perkins  and  Ruedrich  (1973)  and  Sayles  and  Haines  (1974)  and  others. 
Although  the  existence  of  a  true  steady  state  in  the  shorter  duration  tests 
is  questionable.  The  longer  duration  tests  generally  show  creep  exponent 
values  approaching  a  value  of  3.0. 

A  comprehensive  study  of  creep  of  frozen  sands  carried  out  by  Sayles 
(1968)  showed  that  strain  rates  continued  to  attenuate  with  fime  for  tests 


•• 
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carried  out  at  stress  levels  below  the  long  term  strength.  Damped  creep 
behavior  was  also  observed  by  Sayles  and  Haines  (1974)  in  unconfined 
compression  creep  on  fine-grained  remoulded  frozen  samples  when  the  strain 
did  not  exceed  20#  over  the  test  duration. 

The  influence  of  confining  pressure  on  the  creep  rate  of  remoulded 
frozen  Ottawa  sand  was  studied  by  Sayles  (1973).  The  findings  of  this  study 
reported  that  the  amount  of  creep  strain  can  be  reduced  by  increasing  the 
confining  pressure.  Primary  creep  was  still  evident  at  strains  exceeding 
20$.  A  separate  study  carried  out  by  Alkire  and  Andersland  (1973)  showed  an 
exponential  decrease  in  creep  strain  rate  with  an  increase  in  confining 
pressure.  Application  of  confining  pressure  enhances  the  frictional 
behavior  of  the  sand. 

Many  recent  laboratory  studies  have  focussed  attention  on  measuring  the 
empirical  parameters  in  the  constitutive  equation  for  natural  permafrost 
soils.  Thompson  and  Sayles  (1972)  studied  the  in-situ  horizontal  and 
vertical  closure  of  the  rox  Tunnel  in  Alaska.  The  roof  of  the  tunnel  was 
overlain  with  17  m  of  frozen  ice-rich  Fairbanks  silt.  Samples  of  the  silt 
were  tested  in  unconfined  compression  at  stress  levels  between  250  and  2000 
kPa.  The  test  temperature  was  1.7°C.  They  concluded  the  constitutive 
equation  for  the  frozen  silt  has  the  same  form  as  equation  2.4  with  an 
exponent  of  4.  Numerical  analysis  of  the  room  predicted  a  closure  rate  more 
than  3  times  faster  in  the  field  than  that  measured  in  the  laboratory. 

Roggensack  (1977)  tested  ice-rich  fine-grained  g I ac io I acustr i ne  clay 
from  the  Mackenzie  River  Valley.  Triaxial  constant  stress  creep  tests  were 
carried  out  at  deviatoric  stress  levels  between  20  and  400  kPa  and 
temperatures  close  to  -1°C.  This  data  is  shown  in  Figure  2.3.  Although  the 
data  showed  considerable  scatter,  especially  at  stresses  below  100  kPa,  a 
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bilinear  flow  law  was  f  i  t  to  the  data.  A  power  law  fit  to  the  data  above 
100  kPa  gives  an  exponent  of  2.75  which  is  slightly  faster,  but  consistent 
with  the  value  for  polycrystalline  ice.  In  the  low  stress  range,  transient 
processes  play  a  dominant  role  in  the  overall  deformation  response  of  the 
sample.  Increasing  the  confining  pressure  did  not  have  as  a  significant 
effect  on  decreasing  the  creep  rate  as  was  reported  for  coarse-gra i ned  soils 
(e.g.  Sayles,  1973;  Alkine  and  Andersland,  1973). 

McRoberts  et  al,  (1973)  tested  undisturbed  ice-rich  silt  from  Norman 
Wells,  N.W.T.  Constant  stress  creep  tests  were  carried  out  in  a  triaxial 
cel  I  at  deviatoric  stress  levels  between  10  and  690  kPa.  Test  temperatures 
ranged  from  -1  to  -3°C.  Their  results  are  shown  in  Figure  2.4,  adjusted  to 
a  reference  temperature  of  -1°C.  An  upper  bound  to  their  data  is  given  in 
terms  of  a  bi I  inear  flow  law: 


1.6x10-7  1 .5  x  10"14  _ 

•  = - a5  +  -  a6 

( 1 -T ) 1  -  8  d  (1-T)1*8  d 

where  eg  =  axial  strain  rate  (yr-1) 

a ^  =  deviatoric  stress  (kPa) 

T  =  temperature  (°C) 


(2.9) 


A I  I  but  one  of  the  samples  tested  at  a  deviatoric  stress  exceeding  400 
kPa  terminated  in  failure. 
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In  the  low  stress  range  the  value  of  the  creep  exponent  is  3.0, 
consistent  with  polycrystalline  ice.  Tertiary  creep  may  be  influencing  the 
strain  rates  in  the  high  stress  range. 

Savigny  (1980)  tested  fine-grained,  ice-rich  g I ac io I acustr i ne  clay 
samples  from  the  proposed  Arctic  Gas  crossing  of  the  Great  Bear  River, 
N.W.T.  Constant  stress  triaxial  creep  tests  were  carried  out  at  stress  and 
temperature  conditions  simulating  field  behavior.  The  laboratory  results 
reported  showed  no  correlation  with  the  flow  law  for  polycrystalline  ice 
(Morgenstern  et  al.,  1980).  Savigny  (1980)  attributed  this  scatter  to  the 
ground  ice  structure  and  the  generation  of  a  non-uniform  stress  fields 
within  the  specimens  in  response  to  application  of  confining  pressure. 

Nixon  (1978)  reviewed  the  data  of  Roggensack  (1977)  and  McRoberts  et 
al.  (1978).  He  proposed  a  tentative  upper  bound  to  the  creep  rates  of 
McRoberts  et  al.  (1973)  at  a  temperature  of  -2.5°C  and  stress  less  than  100 
kPa: 


c  =  2.135  x  10  8c3 


where  the  units  are  In  (years)  ^  and  kPa. 
2.5  Conclusions 


(2.10) 


The  combined  results  of  Roggensack  (1977),  McRoberts  et  a  I .  (1978)  and 
Savigny  (1980)  are  shown  in  Figure  2.5.  The  two  flow  laws  shown  in  this 
f i gure  are: 


_ 
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1)  flow  laws  for  polycrystalline  ice  (Morgenstern  et  a  I . ,  1980) 

2)  tentative  upper  bound  to  soil  data  of  fine-grained  ice-rich 
soil  proposed  by  Nixon  (1973)  extended  above  100  kPa. 

With  a  few  exceptions,  the  laboratory  data  fall  within  a  narrow  band  at 
stress  levels  exceeding  75  kPa.  Below  this  stress,  transient  processes  may 
still  dom i nate. 

Based  on  the  data  presented  in  Figure  2.5,  it  can  be  concluded  that  the 
creep  law  of  polycrystalline  ice  represents  a  reasonable  upper  bound  to  the 
constitutive  behavior  of  fine-grained,  ice-rich  permafrost.  Exceptions  to 
this  occur  in  the  high  stress  range  where  the  creep  tests  terminated  in 
fa i I ure. 
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FIGURE  2.1 


TYPICAL  CONSTANT  STRESS 
CREEP  CURVES 
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(Bd>l)  SS3y±S  3AI103333 


FIGURE  2.2  COMPARISON  OF  DEFINITION  OF  EFFECTIVE  STRESS 
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FIGURE  2.3  AXIAL  STRESS  VERSUS  STRAIN  RATE 

(Roggensack,  1977) 


FIGURE  2.4  DEVIATORIC  STRESS  VERSUS  STRAIN  RATE 

(McRoberts  etal.,  1978) 
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FIGURE  2.5  STEADY  -  STATE  CREEP  DATA  FOR  ICE  -  RICH  GLACIOLACUSTRINE 

CLAY  FROM  VARIOUS  LOCATIONS  IN  THE  MACKENZIE  RIVER 
VALLEY,  N.W.T.  (Savigny,  1980) 


CHAPTER  I  I  I 


FINITE  ELEMENT  CREEP  ANALYSIS 


3.1 


Genera  I 


The  primary  problem  faced  by  geotechnical  engineers  concerned  with  the 
design  and  construction  of  large  structures  founded  on  permafrost  is  the 
determination  of  the  stress  state  in  the  frozen  soil.  Experimental  evidence 
has  shown  that  time-dependent  deformations  govern  the  behavior  of  frozen 
soils.  The  viscous  nature  of  the  ice  present  in  discrete  layers  and  in  the 
pore  spaces  will  cause  ice-rich  soils  to  creep  at  very  I ow  deviatoric  stress 
I  eve  I s . 

Hoff  (1954)  simplified  the  determination  of  the  stress  state  in  a  non¬ 
linear  viscous  medium  by  proposing  the  elastic  analogue.  The  analogue 
simply  states  that: 

"The  stress  distribution  in  a  body  whose  deformations  are  governed  by  a 
generalized  version  of  a  non-linear  creep  law  is  the  same  as  that  in  a 
non-linear  perfectly  elastic  body  provided  that  the  s tress-stra i n  law 
and  the  boundary  conditions  are  suitably  chosen". 

There  are  two  methods  which  can  be  used  to  solve  creep  problems: 

1)  direct  iteration  procedure  -  the  final  solution  is  obtained  by  an 
iterative  solution  of  the  non-linear  equations 
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2)  incremental  procedure  -  the  final  solution  is  obtained  by  solving 

a  number  of  linearized  problems  between  the  initial  and  final 

time. 

In  the  finite  element  method,  the  inclusion  of  creep  behavior  is  most  easily 
handled  by  using  an  incremental  procedure  and  treating  the  creep  strains  as 
initial  stra i ns . 

The  formulation  of  the  incremental  procedure  was  presented  by  Mendelson 
et  al  (1959).  In  their  approach,  they  extended  the  solution  of  plastic  flow 
problems  using  the  method  of  successive  approximations  to  include  creep 
behavior  with  arbitrary  creep  laws.  In  this  method,  the  elastic  solution  is 

used  to  calculate  the  first  approximation  to  the  creep  strains.  These 

strains  are  then  used  to  calculate  a  new  stress  and  total  strain  state 

within  the  body.  An  improved  creep  strain  approximation  is  then  computed 
and  the  procedure  is  repeated  until  convergence  is  obtained  for  that  time 
i n ter va I.  Essentially,  this  solution  technique  is  i ter at i ve  within  each 
time  increment.  Creep  solutions  were  presented  for  uniaxial  creep  in  a  flat 
plate  and  a  biaxial  case  of  rotating  disks. 

Lin  (1962)  also  presented  an  incremental  procedure  for  analyzing  plates 
subject  to  bending  with  arbitrary  creep  character i sties .  However,  in  this 
approach  the  creep  strain  increments  were  treated  as  equivalent  loads  and 
edge  moments.  This  procedure  is  much  less  complicated  than  the  method 
presented  by  Mendelson  et  al.  (1959)  since  the  method  of  successive 

approximations  may  require  several  interations  within  the  same  time 
i nterva  I  . 

The  application  of  the  finite  element  method  for  solving  non-linear  creep 
problems  using  an  incremental  procedure  where  the  creep  strains  are  treated 
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as  initial  strains  has  been  studied  by  many  researchers  over  the  past  15 

years.  The  incremenTal  initial  strain  procedure  is  widely  used  because: 

it  is  independent  of  the  type  of  creep  law  used; 

2)  it  provides  a  description  of  the  intermediate  stages  of  creep 

it  is  quite  simple  to  extend  basic  finite  element  codes  to  include 

nonlinear  effects  arising  from  creep. 

Gerstenkorn  and  Kobayashi  (1966)  used  the  method  to  solve  for  creep 

deformations  in  a  thick-walled  cylinder  subjected  to  an  internal  pressure. 
Greenbaum  and  Rubinstein  (1968)  developed  a  finite  element  procedure  to 
study  the  time-dependent  behavior  of  axisymmetric  bodies  such  as  thick- 

walled  cylinders  fitted  with  different  types  of  end  closures.  Sutherland 

(1970)  extended  the  work  of  Greenbaum  and  Rubinstein  (1968)  to  include 

problems  subject  to  plane  stress  and  plane  strain  conditions. 


The  incremental  initial  strain  procedure  was  used  by  Emery  (1971)  to 
study  the  creep  behavior  of  slopes  and  excavations  in  rock,  cohesive  soils 
and  ice.  Na  i  r  and  Boresi  (1970)  used  the  incremental  solution  method  to 
study  the  time-dependent  behavior  of  circular  cavities  in  rock  subjected  to 
a  hydrostatic  pressure.  This  procedure  was  also  used  by  Van  Winkel  et  a  I . 

(1972),  to  predict  time-dependent  deformation  of  circular  openings  in  salt 
med i a . 

A  variable  stiffness  method  for  creep  analysis  was  presented  by  Kim  and 
Kuhlemeyer  (1977).  This  method  allows  the  use  of  relatively  large  time 
intervals  for  which  the  initial  strain  procedure  would  become  unstable. 

However,  the  disadvantage  of  this  method  is  the  additional  computer  time 
associated  with  regeneration  of  the  structure  stiffness  matrix  at  the  end  of 
each  time  interval.  This  method  has  proven  useful  in  so i l-structure 
interaction  problems  where  the  structure  exhibits  material  non-linearity. 
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Thompson  and  Sato  (1976)  and  Dawson  and  Tillerson  (1977)  reported  a 
finite  element  formulation  for  simulating  the  creeping  flow  of  an 
incompressible  material.  The  stress-stra i n  rate  response  of  the  material 
can  either  be  linear  (i.e.  Newtonian  fluid)  or  non-linear  (i.e.  non- 
Newtonian  fluid).  This  method  allows  the  ability  to  follow  the  flow  through 

very  large  changes  in  material  geometry.  Velocity  components  are  chosen  as 
the  primary  variables. 

The  velocity  field  is  constrained  by  the  material  incompressibility 
condition.  This  constraint  is  incorporated  into  the  potential  energy 
functional  through  the  use  of  a  Lagrange  multiplier.  The  multiplier  is 
physically  interpreted  as  the  mean  normal  stress  or  the  pressure  within  the 
element.  Thus,  approximating  functions  for  the  pressure  as  well  as  the 
velocity  for  each  element  must  be  chosen.  This  adds  greatly  to  the  solution 
time  required  to  solve  the  field  equations. 

Numerical  difficulties  arise  when  trying  to  satisfy  the  material 
incompressibility  constraint.  To  ensure  complete  incompressibility 
everywhere  within  the  element,  the  approximating  functions  chosen  for  the 
pressure  within  the  element  must  be  the  same  degree  as  that  used  for  the 
element  strain.  Thus,  if  a  linear  approximation  to  the  pressure  field  is 
chosen,  a  full  quadratic  approximating  function  must  be  chosen  for  the 
components  of  the  velocity  field.  The  complete  incompressibility 
restriction  on  the  velocity  field  has  been  shown  to  be  give  slow  convergence 
for  coarse  element  grids. 

Material  incompressibility  can  be  satisfied  in  an  average  sense  over 
the  entire  element.  This  is  accomplished  by  assuming  a  constant  pressure 
and  linear  variation  for  the  strain  rates  and  dilatation  within  each 
element.  In  this  approach,  the  total  volume  of  the  element  remains 
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unchanged  although  local  incompressibilities  within  the  element  are  still 
free  to  exist.  Relaxation  of  the  incompressibility  constraint  will  give 
more  accurate  results  for  coarse  element  layouts.  However,  the  constant 
pressure  assumption  requires  a  fine  element  layout  in  situations  where  an 
accurate  evaluation  of  element  stresses  are  desired. 

3*2  Incremental  Initial  Strain  Procedure 


As  a  material  undergoes  creep,  the  stresses,  a  j  j »  in  the  body 
will  change  with  time  as  shown  schematically  in  Figure  3.1.  In  the 
incremental  procedure,  the  smooth  stress-time  curve  is  replaced  by  a  series 
of  incremental  steps.  Each  step  consists  of  a  constant  stress  period, 
At,  followed  by  an  instantaneous  increment  of  stress. 

The  fundamental  assumption  in  the  incremental  creep  solution  technique 
is  that  the  time  can  be  subdivided  into  sufficiently  smal  I  time  intervals 
such  that  the  stress  can  be  assumed  constant  within  each  time  interval.  If 
valid,  the  non-linear  creep  problem  can  be  solved  as  a  series  of  linear 
problems  for  each  time  interval  using  the  stresses  of  the  previous  time 
interval  to  calculate  new  increments  of  creep  strain  and  treating  them  as 
initial  strains  for  the  current  time  period. 

The  following  assumptions  are  made  in  order  to  establish  a  stress- 
strain  relationship  for  a  creeping  material: 

1)  there  is  no  volume  change  associated  with  the  creep  strains 

2)  for  an  isotropic  material  the  principal  directions  of  the  strain 
rate  and  stress  tensors  coincide; 

the  creep  rate  is  independent  of  any  superimposed  hydrostatic 
state  of  stress; 
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4)  the  generalization  of  the  uniaxial  creep  laws  to  the  multiaxial 
state  of  stress  should  recover  the  uniaxial  relationship  for  the 
case  of  uniaxial  stress. 

Assumptions  1,  2  and  3  can  be  stated  mathematically  as: 


(3.1  ) 


a  S  .  . 

'  J 


where:  a  =  constant  of  proportionality 

This  equation  states  that  the  increments  of  creep  strain  rate  are 
proportional  to  the  i ns i antaneous  values  of  the  deviatoric  stress  tensor  and 
are  independent  of  stress  history. 

Substituting  equation  3.1  into  Odqvist's  definition  of  effective  strain 
rate  (i.e.  equation  2.3)  and  making  use  of  the  definition  of  effective 
stress  yields: 


a  =  3 (Ae°/a  ) 


(3.2) 


e  e 


Substituting  equation  3.2  into  3.1  gives  the  stress-stra i n  rate  relationship 
for  a  creeping  material: 


(3.3) 


. 
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For  the  plane  strain  case,  the  following  equations  are  valid: 


2  2  2  ?  ? 
cr  "a  )  +  ( a  “a  )  +  (a  “a  )  ”  +  6x 

e  2  x  V  Y  z 


z  x 


xy 


(3.4) 


til  =  4[<aJC>2  +  (A*C)2  +  4|CAJC  +  ( ay  /2)2 
3  x  y  x  y  xy 


(3.5) 


C  C 

Ae  =  (Ae  /a  ) (2a  -a  _a  )/2 
x  e  e  x  y  z 


(3.6) 


C  C 

Ae  =  (Ae  /a  ) (2a  -a  ~a  )/2 
y  e  e  y  z  x 


(3.7) 


C  C 

Ae  =  (Ae  /a  )  (2a  ~a  ~a  )  /2 
z  e  e  z  x  y 


(3.8) 


AyC  =  3( Ae°/a  )T 
xy  e  e  xy 


(3.9) 


The  essential  feature  of  the  incremental  procedure  is  to  proceed  in 
small  intervals  of  time  and  relate  an  increment  of  strain  to  an  increment  of 
stress.  The  final  stress  and  strain  state  in  the  body  is  obtained  by 
summing  each  increment. 

At  the  beginning  of  each  time  interval,  the  stresses,  elastic  strains 
and  creep  strains  are  known  from  the  calculations  of  the  previous  time 
increments.  At  time  t=0,  the  elastic  solution  is  used  as  the  starting 
point.  The  incremental  procedure  for  each  increment  of  time  can  be 
summarized  by  the  following  five  steps: 
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1)  obtain  the  value  of  the  effective  stress  (ae)  from  the  stresses 
(aj;)  determined  from  the  previous  time  step; 

2)  calculate  the  effective  creep  strain  increment  (Ae  )  by  substi- 

©  ' 

tuting  the  effective  stress  (a0)  (determined  from  step  1)  into  the 
appropriate  creep  relationship; 

3)  calculate  creep  strain  increments  (Aejj)  in  each  direction  with 
the  stresses  ( cr  j  j )  of  the  preceding  time  step,  the  effective 

stress  (a^)  (from  step  1)  and  the  effective  creep  strain  increment 
^se)  (determined  in  step  2); 

q 

4)  treating  the  creep  strain  increments  (Aejj)  (calculated  in  step  3) 

as  initial  strains  and  using  the  constitutive  equations,  boundary 
conditions  and  equilibrium  equations  for  the  particular  problem, 
calculate  the  increment  of  stress  (  Aa^j  )  at  the  end  of  the  time 
interval;  (the  initial  strains  are  converted  to  ficticious  creep 
forces  applied  at  the  nodal  points); 

y)  the  incremental  stresses  (Ae.  .)  (obtained  in  step  4)  are  added  to 

the  stress  of  the  previous  time  increment  to  obtain  a  new  stress 
distribution  for  the  current  time  interval. 

Tnis  procedure  is  repeated  for  each  time  interval  until  either  the 
final  time  is  reached  or  until  the  stress  distribution  does  not  change  (i.e. 
a  steady  state  condition  has  been  achieved).  As  long  as  the  incremental 

stresses  are  small  compared  to  the  previous  stresses,  the  basic  assumption 
of  i ne  incremental  procedure  is  not  violated.  General ly,  this  condition  can 
be  met  by  selecting  time  increments  small  enough  to  yield  the  desired 
accuracy  without  employing  iterative  techniques. 
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3.3 


Cumulative  Creep  Rule 


In  general,  the  effective  incremental  creep  strain  rate  is  a  function 
of  the  effective  stress,  total  effective  creep  strain,  temperature  and  the 
strain  history  of  the  material.  Thus,  a  general  functional  form  of 
is  expressed  as  follows: 


(3.10) 


where:  T  =  temperature 

t  =  time 


Materials  exhibiting  damped  creep  response  have  a  strain  rate  that  is  a 
function  of  time  as  well  as  applied  stress.  An  expression  must  be  obtained 
from  constant  stress  creep  tests  that  is  able  to  predict  the  strain  rate 
under  varying  states  of  stress.  Assume  the  damped  constant  stress  creep 
data  can  be  presented  by  an  empirical  relationship  of  the  form: 


(3.11 ) 


where:  m  =  time  exponent 


The  creep  strain  rate  can  be  obtained  by  differentiating  equation  3.11 
with  respect  to  time: 


e 


(3.12) 
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where:  f(a  t) 


general  function  form  expressing  the  creep  strain  rate 
as  a  function  of  applied  stress  and  time. 


The  time  variable  can  be  el iminated  from  the  creep  strain  rate 
expression  by  substituting  equation  3.11  into  equation  3.12: 


e 

e 


mA 1/m  n/m  (1-1/m) 
mA  o  e 
e  e 


(3.13) 


general  function  form  expressing  the  creep  strain 
rate  as  a  function  of  applied  stress  and  total  creep 
stra i n . 

It  is  clear  that  integration  of  equation  3.12  and  3.13  will  not  give 
the  same  result  for  the  general  case.  Two  basic  rules  that  are  currently 
used  to  describe  the  time  dependency  of  strain  for  creep  in  metals,  plastics 
and  concrete  are  the  strain  hardening  law  and  the  time  hardening  law. 


where:  g(ae,ee)  = 


The  strain  hardening  law  for  a  set  of  isothermal  constant  stress  curves 
is  illustrated  in  Figure  3.2.  This  hypothesis  implies  the  creep  strain  rate 
is  a  function  only  of  the  instantaneous  value  of  stress  and  the  total 
accumulated  creep  strain.  A  change  in  stress  state  is  represented  by  a 
horizontal  line  from  one  constant  stress  curve  to  the  new  stress  curve  at 
the  same  total  creep  strain.  The  total  time  is  given  by  the  sum  of  al  I 
constant  stress  time  intervals.  This  concept  assumes  a  mechan i ca I -equat i on- 
of-state  exists  for  a  creeping  material  relating  the  variables  stress, 
strain  and  strain  rate.  The  behavior  is  independent  of  the  stress  history 
the  material  was  subjected  to  in  the  early  stages  of  creep.  This  is 
analogous  to  the  theorem  proposed  by  Odquist  (1966)  which  stated: 
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"If  a  test  piece  is  subjected  to  a  series  of  positive  stress  values, 
a1  *  °2>  •••>  an»  each  acting  during  a  time  period  At^  ,  At£, 

^n*  ^hen  The  resulting  total  creep  is  independent  of  the  order  in 
which  ; he  stress  values  were  app! ied,  each  during  its  respective  time 
per iod" . 


The  time  hardening  concept  for  a  set  of  isothermal  constant  stress 
creep  curves  is  illustrated  in  Figure  3.3.  This  cumulative  rule  assumes  the 
creep  strain  rate  is  governed  by  the  instantaneous  value  of  the  stress  and 
total  time  from  the  beginning  of  the  test  and  is  independent  of  the  stress 
history  of  the  material  prior  to  any  point  in  time.  A  change  in  stress 
stafe  is  represented  by  a  vertical  I  ine  from  one  constant  stress  curve  to 
the  new  stress  curve  at  the  same  total  elapsed  time.  The  total  accumulated 

creep  strain  is  given  by  the  sum  of  the  creep  strain  increments  during  each 
time  i nterva I . 

The  validity  of  these  two  hypothesis  have  not  been  specifically  studied 
in  creep  tests  on  permafrost  samples.  Experience  gained  in  metal  and 
plastic  creep  as  well  as  undrained  creep  behavior  of  cohesive  soils  suggests 
that  the  strain  hardening  cumulative  creep  law  adequately  describes  the 
creep  strain  rate  under  varying  stress  conditions.  Ladanyi  and  Johnston 
(1973)  used  the  strain  hardening  hypothesis  to  determine  the  in-situ  creep 
parameters  for  a  frozen  varved  silty  clay  from  multistage  pressuremeter 
data.  The  strain  hardening  cumulative  creep  law  can  also  be  used  to  predict 
the  creep  rate  of  a  material  under  varying  temperature  conditions  (Dorn, 
1961 ) . 

Derivation  of  a  cumulative  creep  law  which  included  all  the  factors 
governing  the  creep  rate  of  a  material  would  be  a  formidable  task. 
Cumulative  creep  laws  have  been  put  forward  in  the  metal  creep  literature 
which  take  stress  and  temperature  history  into  account.  However,  these  can 
lead  to  untractanle  mathematical  expressions. 
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—in i t8  Element  Formulation  of  Incremental  Initial  Strain  Method 


The  basis  of  the  finite  element  method  is  to  represent  the  continuum  as 
an  assemblage  of  finite  elements.  These  elements  are  interconnected  at  a 
finite  number  of  points  referred  to  as  nodal  points.  Approximating 
functions  are  chosen  to  represent  the  variation  of  the  field  variable  over 
the  element.  A  variational  principle  of  mechanics  is  used  to  derive  the 
equilibrium  equations  for  each  element.  The  equilibrium  equations  for  the 
entire  body  are  obtained  by  summing  all  the  individual  element  equations 
with  proper  regard  to  displacement  continuity  at  the  nodal  points.  The 
boundary  conditions  are  applied  and  the  entire  set  of  equations  are  solved 
to  obtain  the  unknown  displacements  at  the  nodal  points. 


in  i he  present  study,  the  Theorem  of  Minimum  Potential  Energy  is  used 
to  formulate  the  equilibrium  equations.  The  elements  used  are  constant 
strain  triangles.  Within  each  element,  displacement  functions  are  chosen 
such  that  inter-element  compatab i I i ty  is  maintained.  The  displacements  are 
assumed  to  be  linear  functions  of  the  coordinates.  The  potential  energy  of 
an  element  is  calculated  by  subtracting  the  work  done  by  the  external  forces 
from  the  stored  strain  energy.  Minimization  of  the  potential  energy  yields 
tne  equilibrium  equations  for  the  element.  The  equilibrium  equations  for 
i he  entire  assemblage  are  obtained  by  summing  the  individual  element 
contributions.  The  resulting  set  of  linear  algebraic  equations  are  solved 
for  the  unknown  displacements.  A  detailed  derivation  of  the  equilibrium 
equations  for  the  constant  strain  triangle  element  can  be  found  in  numerous 
standard  finite  element  texts  [see  e.g.  Desai  and  Abel  (1972)  and 
Zienkewicz  ( 1977) ] . 


The  creep  problem  is  solved  by  assuming  the  change  in  total  strain 
during  one  time  interval  is  composed  of  a  change  in  elastic  and  creep 
stra ins. 


. 
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TEC  r 

Ae  -  Ae~  +  Ae  =  I [ Aa  -v(A a  +A a  )]  +  Ae° 
Y  x  x^y  2  x  y 


(3.14) 


TEC 

Ae  -  Ae  +  Ae  -  i [ Aa  -v(Aa  +Aa  ) ]  +  Ae 
z  x  x  T  2  x  y 


C 


z 


where  the  superscript  T,  E  and  C  denote  total,  elastic  and  creep 
respect i ve I y. 

Plane  strain  problems  have  a  kinematic  constraint  placed  on  the  strain 
in  Tie  Z  direction  (i.e.  Ae^"z=o).  Thus,  the  third  expression  in  equation 
3.14  gives  the  value  for  the  stress  in  the  Z  direction: 


C 


v(Aa  +Aa  )  -  EAe 

x  y 


(3.15) 


z 


z 


Substituting  equation  3.15  into  the  first  two  expressions  in  equation  3.14 
yields: 
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I  _  i  C  C 

Ae  -  I  [  ( 1  —  v  )Aa  -v(1+v)Acr  1  +  Ae  +  vAe 
x  P  x  y  x  z 


(3.16) 


?  r  r 

Ae  =  It  ( 1 -v  )Act  -v(1+v)Aa  1  +  As  +  vAe 

y  E  y  x  y  z 


Equation  3.16  shows  that  to  account  for  plane  strain,  a  strain  equal  to 
q 

vAez  must  be  added  to  the  strains  in  the  X  and  Y  directions.  The 
increments  of  creep  strain  Ac?.  arg  ca|cu|afed  from  equa+ion  3<3- 

The  strain-displacement  relations  and  application  of  the  Theorem  of 
Minimum  Potential  Energy  yields  the  equilibrium  equations  for  one  time 
interval  as: 


[ K ] { Aq}  =  {p}  +  (F  } 

c 


(3.17) 


where:  IK]  =  structure  stiffness  matrix 

{Aq}  =  incremental  nodal  displacements 
{P}  =  nodal  load  vector 

<fy  =  vector  of  ’fictitious*  nodal  loads  due  to  an  increment  of 
creep  strain. 


The  introduction  of  creep  strain  increments  into  the  equilibrium  equations 
using  the  initial  strain  method  is  given  in  Appendix  A. 
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At  the  end  of  each  time  interval,  the  vectors  { P}  and  {F  }  and 
boundary  conditions  are  known.  Equation  3.17  is  solved  for  the  unknown 
displacement  increments  by  the  Choleski  square  root  method  for  symmetrically 
banded  matrices.  For  linear  elastic  analysis,  the  structure  stiffness 
matrix  does  not  change  with  time  and  need  only  be  generated  once  for  the 
initial  time  increment  and  stored  to  be  used  for  each  successive  time 
increment.  Should  the  magnitude  of  the  displacements  become  excessive,  the 
geometry  may  need  to  be  updated  and  the  stiffness  matrix  regenerated  to 
reflect  the  new  grid  geometry.  This  was  not  necessary  for  the  small 
displacement  problems  considered  in  this  research. 

3.5  Time  Increment  Selection 


The  selection  of  a  time  increment  to  ensure  the  solution  process  remain 
stable  is  very  important  in  the  incremental  procedure.  During  the  early  or 
transient  stages  of  creep,  the  stresses  are  changing  very  rapidly  so  it  is 
necessary  to  select  sma  I  I  time  increments.  As  the  solution  approaches  a 
steady  state  condition,  the  magnitude  of  the  time  interval  can  be  increased. 

Experience  of  Greenbaum  and  Rubinstein  (1968),  Sutherland  (1970)  and 
Emery  (1971)  shows  that  the  solution  process  becomes  unstable  if  the  maximum 
change  in  creep  strain  is  larger  than  the  effective  creep  strain.  The 
magnitude  of  time  intervals  chosen  is  a  function  of  the  stress  state  in  the 
body  and  the  specific  form  of  the  creep  law. 

The  initial  time  increment  is  chosen  such  that  the  maximum  ratio  of  the 
effective  creep  strain  increment  to  the  effective  elastic  strain  is  equal  to 
or  I  ess  than  i .e. : 


. 
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r  C  ,  E, 

[Ae  /Ae  ]  <  n 

e  o  max  o 


(3.18) 


where:  .04  <  n0  <  .10 


Successive  time  increments  are  computed  on  the  basis  of  the  fractional 
change  in  maximum  effective  stress  allowed  per  time  interval,  i.e.: 


[ w/( Act  /a  )  ]  At . 

e  e  max  i 


(3.19) 


where:  i  )  .035  <  to  <  .10 

or  ii)  1.2  <  u/(Act  /or  )  <2.0 

e  e  max 

The  maximum  ratio  of  effective  creep  strain  increment  to  the  effective 
elastic  strain  for  successive  time  increments  must  also  satisfy: 


[Ae^/Ae^]  < 

e  e  max 


(3.20) 


where:  0.1  <  n-j  <  1»0 

The  optimum  values  for  nQ,  n  ^  and  u  are  a  function  of  the  specific 
creep  law  and  the  stress  state  that  exists  in  the  body.  Thus,  it  is 
advantageous  to  compute  the  value  for  the  time  interval  internally. 

The  manner  in  which  the  program  calculates  the  time  interval  is  as 
fo I  I ows : 


Assume  for  generality  the  creep  behavior  of  the  material  is 
governed  by  the  creep  relationship  given  by  equation  3.11,  i.e.: 
(the  steady  state  creep  law  is  derived  from  this  expression  by 
sett i ng  m  =  1.0),  i.e.: 


1) 
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The  initial  time  interval  is  computed  by  combining  equation  3.18 
with  the  above  primary  creep  relationship: 


1  /  m 


/AE I 


Ae 


(3.21 ) 


where:  E  =  Young's  modulus 

The  effective  stress  is  calculated  based  on  the  elastic  stress 
distribution  at  t  =  0. 

Succeeding  time  intervals  are  obtained  by  computing  the  maximum 

ratio  (  e/ae)  and  satisfying  the  conditions  given  by  equation 
3.19.  A  further  restriction  on  succeeding  time  intervals  to 
prohibit  divergence  is  satisfaction  of  equation  3.20. 


^  ficticious  time,  t^,  is  introduced  for  materials  governed  by  a 
primary  creep  law  and  obeying  a  strain  hardening  cumulative  creep  law.  The 
fictitious  time  is  required  to  calculate  the  effective  creep  strain 

increments  for  all  increments  i  (where  i  =2,  3,  . ..,  n)  total  number  of 
i ncrements) : 


Ae  ^  =  Aan I ( t ,+a t) m 
e  e  f 


(3.22) 


where  the  fictitious  time  is  given  by: 


(3.23) 
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The  ficticious  time  concept  is  illustrated  in  Figures  3.4  for  an 
increasing  and  decreasing  stress  case. 


Once  a  steady  state  condition  is  reached  the  stresses  do  not  change 
with  time  and  the  creep  strains  can  be  extrapolated  to  the  final  time.  Two 
tests  can  be  made  to  check  for  the  existence  of  steady  state  creep: 

1)  Check  the  maximum  fractional  change  in  effective  stress,  i .e. : 


t A a  /a  1  < 

e  e  max 


(3.24) 


2)  Check  the  maximum  change  in  effective  stress  per  unit  time,  i.e.: 


t  Act  /a "I"  •  1  <  X 

e  i  max 


(3.25) 


Values  of  p  and  x  depend  on  the  specific  problem  to  be  solved. 

The  steady  state  condition  checks  can  also  be  used  to  monitor  the 
stability  of  the  solution  procedure.  Erratic  fluctuations  of  the  maximum 
fractional  change  in  effective  stress  or  the  maximum  change  in  effective 
stress  per  unit  time  would  indicate  solution  instability.  To  alleviate  this 
problem,  smaller  values  of  X]^t  ^  ^ ,  and  cj  are  required. 

3.6  Finite  Element  Programme 

The  finite  element  programme  for  time  dependent  creep  deformation 
analysis  was  originally  developed  by  Emery  (1971)  at  the  University  of 
British  Columbia.  Minor  changes  were  made  to  the  programme  after  it  was 
obtained  by  the  University  of  Alberta. 


* 
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The  incremental  initial  strain  method  described  in  a  previous  section, 
is  used  to  solve  the  creep  response.  The  programme  is  capable  of  solving 
isothermal  creep  problems  with  linear  elastic  material  behavior  in  plane 
strain  situations  only.  The  element  type  used  is  the  constant  strain 
triangle.  External  loads  consist  of  nodal  point  loads  and/or  gravity 
I oad i ng. 

The  general  form  of  the  flow  law  implemented  into  the  programme  is  given  by: 


nl  n? 

[Ve  +  Ve  1+m 


(3.26) 


where: 


ee  -  effective  strain  rate 
ae  =  effective  stress 
t  =  time 


A-],A2  =  creep  coefficients 

n i , n 2  =  creep  stress  exponents 

m  =  creep  time  exponent 


The  parameters  ,  n 2  and  m  are  material  constants.  This 

general  form  of  flow  law  was  selected  because  it  degenerates  to  a  simple 
power  law,  a  bilinear  flow  law  and  a  primary  (i.e.  transient)  power  law. 


A  stop  and  restart  capability  has  been  implemented  into  the  finite 
element  program.  This  feature  allows  the  creep  solution  to  be  interrupted 
after  a  specified  time  or  number  of  increments.  Solution  results  up  to  that 
point  may  be  examined,  or  material  properties  and/or  externally  applied 
loads  altered  to  reflect  changes  in  the  creep  simulation. 


A  users'  manual  for  the  program  CREEP  is  presented  in  Appendix  C.  The 
programme  listing  is  given  in  Appendix  D.  A  comparison  of  the  closed-form 
thick  wall  cylinder  solution  with  the  finite  element  results  is  presented  in 
Appendix  E. 
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TIME  t 


FIGURE  3.1 


APPROXIMATION  OF  A  SMOOTH  STRESS-TIME 
CURVE  FOR  INCREMENTAL  PROCEDURE 


STRAIN  STRESS 
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A 


FIGURE  3.2 


STRAIN  HARDENING  CUMULATIVE  CREEP  RULE 


STRAIN  STRESS 
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TIME 


FIGURE  3.3 


TIME  HARDENING  CUMULATIVE  CREEP  RULE 


CREEP  STRAIN  CREEP  STRAIN 
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A 


FIRST  TIME  INTERVAL 


a)  INCREASING  STRESS 


A 


b)  DECREASING  STRESS 


FIGURE  3.4 


CALCULATION  OF  INCREMENT  OF  CREEP  STRAIN 
FOR  CHANGING  STRESS  CONDITIONS 


‘ 


CHAPTER  IV 


IN  SITU  CREEP  ANALYSIS  OF  A  SLOPE  IN  ICE-RICH  PERMAFROST 


4.1  General 


The  proposed  Arctic  Gas  pipeline  route  crosses  the  Great  Bear  River  7.2  km 
above  the  confluence  of  the  Mackenzie  and  Great  Bear  Rivers  (see  Figure 
4.1).  The  left  bank  of  the  Great  Bear  River  is  consistently  steeper  than 
the  right  bank  for  several  kilometres  upstream  or  downstream  from  the 
crossing.  Naturally  occurring  creep  was  a  concern  at  this  particular  site 
because  the  left  bank  rises  at  an  average  slope  22  degrees  over  46  m  from 
river  level  to  the  local  plain  level  (see  Figure  4.2).  Mean  river  elevation 
is  66  m  a.m.s.l.  The  local  plain  elevation  is  102  m  a.m.s.l.  The  top  of 
the  plain  has  very  little  topographic  relief  south  of  the  slope  crest. 

The  left  bank  of  the  Great  Bear  River  at  the  proposed  Arctic  Gas  crossing 
was  studied  by  Savigny  (1980).  An  extensive  site  investigation  programme 
including  drilling,  sampling  and  instrumentation  was  carried  out. 
Undisturbed  permafrost  core  samples  were  obtained  for  high  quality 
laboratory  testing.  The  instrumentation  was  carefully  monitored  for  a  two 
year  period.  This  chapter  presents  the  results  of  a  finite  element 
prediction  of  the  in-situ  creep  behavior  of  the  proposed  Arctic  Gas  crossing 
of  the  Great  Bear  River. 
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4.2 


Site  Geology 


The  sequence  of  geologic  sediments  encountered  at  the  boreholes  drilled 
at  tnis  site  are  shown  in  Figure  4.3.  The  slope  is  composed  mainly  of  sand 
and  clay.  Till  outcrops  near  river  level.  The  river  has  formed  in 
mterbedded  clay  and  sand  which  overlies  shale  and  siltstone  bedrock. 

The  shale  and  siltstone  bedrock  are  highly  weathered  and  soft.  These 
rocks  can  be  likened  to  an  overconsolidated  soil.  The  bedrock  is  overlain 
by  an  alluvial  deposit  of  interbedded  silty  clay  and  sand  approximately 
6.7  rn  thick.  The  two  major  components  of  this  unit  are  highly  plastic  silty 

clay  and  clayey  silt.  Ground  ice  is  present  in  either  reticulate  or 
segregated  forms. 

A  dense  glacial  till  of  Wisconsin  age  lies  uncon formab I y  on  the 
alluvial  deposits.  The  soil  matrix  consists  of  a  low  to  medium  plastic 
sand-silt-clay  mixture.  The  till  is  highly  fissured  where  it  outcrops  along 

tic  va  I  I  oy  ,/a  i  I  .  Pore  ice  and  reticulate  ice  are  the  two  most  common  types 
of  ground  ice  present. 

A  fine-grained  g I ac iol acustr ine  clay  rests  on  the  till  and  is  overlain 
by  g lac iodel ta ic  sand.  The  g I ac io I acustr i ne  clay  was  deposited  in  ice 
dammed  lakes  in  late  Wisconsin  time.  The  soil  matrix  consists  of  a  medium 
to  highly  plastic  silty  clay.  This  deposit  is  approximately  18  m  thick  at 
the  proposed  crossing  site.  The  g I ac io lacustr ? ne  clay  is  ice-rich. 
Reticulate  ice  is  most  common  with  primary  vertical  veins  and  secondary 
horizontal  veins.  Ice  veins  up  to  20  cm  thick  have  been  observed. 

Stratified  ice  is  also  common  at  the  till  contact.  The  ice  content  of  this 
unit  increases  where  the  clay  outcrops  along  the  valley  wall. 


. 

'%  *,  V  tr.  £  *9  V 

' 


i  ; 


% 


50 


The  uppermost  unit  at  this  site  ?  s  a  20  m  thick  stratum  of 
g I ac i ode  I ta i c  sand.  This  material  was  deposited  fol lowing  the  disappearance 
of  the  large  ice  dams  upstream  on  the  Mackenzie  River  near  Fort  Good  Hope, 
N.W.T.  The  sand  is  medium  to  fine  grained  with  thin  layers  of  low  plastic 
silt  common  throughout.  Excess  pore  ice  is  the  dominant  type  of  ground  ice. 

Some  steeply  dipping  ice  veins  are  also  present  but  not  to  any  great 
extent.  ihe  thickness  of  the  active  layer  in  the  sand  varies  anywhere  from 
.3  to  5  m. 


4.3  Field  Instrumentation 

The  field  instrumentation  installed  on  the  left  bank  of  the  Great  Bear 
River  at  the  proposed  pipeline  crossing  consisted  of: 

1)  borehole  inclinometers  to  record  in-situ  downslope  creep 
velocities  in  ice-rich  permafrost  soils  (boreholes  GB1 A,  GB2  and 
GB3 ) ; 

2)  thermistors  to  measure  the  in-situ  geothermal  gradient  (boreholes 
GB1 A,  G32  and  GB3 ) ; 

3)  piezometers  to  measure  pore  pressures  below  the  base  of  the 
permafrost  (boreholes  GB1B  and  GB3A). 

The  location  of  all  the  boreholes  is  shown  on  Figure  4.2. 

The  pore  pressures  measured  at  borehole  GB3A  show  a  close 
correspondence  between  river  level  elevation  and  hydrostatic  pore  pressure. 
This  is  to  be  expected  because  the  coarse  alluvium  that  lies  on  top  of  the 
bedrock  at  this  site  provides  a  free  drainage  path. 
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The  ground  temperature  contours  below  the  depth  of  zero  mean 
temperature  fluctuation  for  the  left  bank  of  the  Great  Bear  River  at  the 
proposed  pipeline  crossing  are  shown  in  Figure  4.4.  These  contours  were 
inferred  from  the  thermistors  installed  in  boreholes  GB1A,  GB2  and  GB3.  The 
temperature  data  indicate  that  the  permafrost  at  this  site  is  warm  (i.e. 
warmer  than  -3°C) . 

At  the  top  of  the  slope  (GB1A)  the  active  layer  is  approximately  3  m 
thick  and  the  depth  of  annual  zero  mean  temperature  fluctuation  is  9  to 
10  m.  The  depth  of  permafrost  inferred  from  the  geothermal  gradient  is  at 
least  61  m.  The  temperature  readings  over  the  two  year  period  show  that 
removal  of  the  surface  vegetal  cover  has  initiated  a  warming  trend. 

The  thermal  regime  of  the  slope  differs  from  that  on  top  of  the  slope. 
This  is  primarily  due  to  the  aspect  of  the  slope  and  the  difference  in 
vegetation.  Along  the  slope  the  active  layer  is  approximately  1.2  m  thick 
and  the  depth  of  annual  zero  mean  temperature  fluctuation  is  6  to  7  m.  At 
approximately  mid-height  of  the  slope,  the  mean  ground  surface  temperature 
•s  -3.3°C  and  near  the  toe  of  the  slope,  the  warming  effect  of  the  Great 
Sear  River  increases  the  mean  ground  surface  temperature  to  -2.5°C.  The 
steeper  thermal  gradient  measured  near  the  toe  of  the  slope  (borehole  GB3) 
is  also  a  result  of  the  thermal  disturbance  offered  by  the  Great  Bear  River. 
The  base  of  the  permafrost  rises  to  the  river  bed  elevation  as  permafrost  is 
absent  below  the  Great  Bear  River. 

Three  inclinometers  were  installed  to  monitor  downslope  creep 
movements : 


.* 
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i)  incl  inometer  GB1A  was  ins+al  led  at  the  top  of  the  slope  through 
the  g I ac i of  I uv i a  I  sand  and  g I ac io I acustr i ne  clay  and  terminates  in 
the  silty  clay  till; 

ii)  inclinometer  GB2  was  installed  at  mid-height  of  the  slope  where 
the  g I ac io I acustr i ne  clay  outcrops  along  the  valley  wall  and 
terminates  at  the  base  of  the  till; 

iii)  inclinometer  GB3  was  installed  near  the  toe  of  the  slope  through 
the  g  I  ac  io  I  acustr  i  ne  clay  and  till  and  terminates  at  the  base  of 
the  interbedded  clay  and  sand. 

The  location  of  each  inclinometer  is  shown  in  Figure  4.5.  The 
inclinometers  were  read  on  twelve  occasions  over  a  two  year  period  extending 
from  May  6,  1975  to  June  14,  1977. 

The  inclinometer  data  show  a  strong  correlation  between  deformations 
and  ice  lenses  in  the  g I ac io I acustr i ne  clay.  In  zones  containing  large  ice 
lenses  separated  1  to  2  m  apart,  the  movements  were  large  and  abrupt  and  are 
concentrated  at  the  ice  lense  location.  In  zones  where  the  ice  lenses  are 
spaced  closer  than  1  m,  the  movements  are  generally  smaller  and  more  gradual 
through  the  zone.  No  net  transverse  slope  movements  of  any  practical  value 
were  recorded. 

Movements  recorded  in  the  overlying  sand  are  much  more  inconsistent 
than  those  measured  in  the  clay.  Savigny  (1980)  attributes  this  to  drilling 
difficulties  encountered  in  the  granular  material  causing  disturbed  zones 
around  the  hole  which  lead  to  non-uniform  stress  distributions  in  the  frozen 
sand. 
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The  velocity  versus  depth  profiles  for  all  three  inclinometer  locations 
are  shown  in  Figure  4.5.  A  net  downslope  velocity  was  recorded  through  the 
medium  to  high  plastic  g I ac i o I acustr i ne  clay  at  hole  GB1 A .  The  maximum 
velocity  was  recorded  near  the  sand-clay  contact.  Erratic  movements  were 
recorded  between  the  29  m  and  34  m  depth. 

Examination  of  the  GB1 A  borehole  stratigraphy  reveals  that  this  depth 
interval  corresponds  to  a  zone  containing  large  pervasive  ice  lenses  spaced 
greater  than  1  m  apart.  The  velocity  maxima  within  this  zone  coincide  with 
the  location  of  large  ice  lenses.  Above  the  29  m  depth  and  below  the  34  m 
depth,  there  is  a  pseudo- I i near  increase  in  velocity  with  increasing 
distance  away  from  the  base  of  the  clay  layer.  The  steeper  velocity 
gradient  above  29  m  i s  a  result  of  large  ice  lenses  spaced  closely  together. 
Below  34  m,  the  ice  lenses  are  much  smaller. 

A  uniform  downslope  velocity  was  recorded  in  the  upper  12  m  of  the 
overlying  sand  at  hole  GB1A.  The  magnitude  of  this  velocity  is 
approximately  one  half  of  the  velocity  recorded  in  the  clay  just  below  the 
sand/clay  contact.  The  12  to  20  m  depth  interval  of  the  sand  predominantly 
shows  an  upslope  velocity.  The  abrupt  change  at  12  m  from  upslope  to 
downslope  velocity  could  possibly  be  due  to  a  large  ice  vein.  At  borehole 
GB1 ,  a  large  ice  vein  was  observed  13  m  below  the  surface. 

No  net  creep  movements  were  recorded  in  the  transverse-to-s I  ope 
direction  at  the  GB1A  inclinometer. 

The  velocities  recorded  through  the  medium  to  high  plastic  clay  are 
very  erratic  at  the  GB2  inclinometer  location.  Savigny  (1980)  describes  two 
causes  for  the  erratic  in-situ  velocity: 
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1)  difficulties  encountered  in  installation  and  monitoring  of 
i nstrumentation 

2)  existence  of  in-situ  sheared  material  as  evidenced  by  the  numerous 
slickensides  in  the  core  samples  retrieved  from  this  hole. 

A  small  net  downslope  creep  movement  is  evident  in  the  silty  clay  till 
at  the  G32  location.  There  is  a  pseudo- I i near  increase  in  creep  velocity 
from  the  base  of  the  til!  up  to  the  t i I  I /g I ac io I acustr i ne  clay  contact. 

A  net  downstream  creep  velocity  was  recorded  in  the  g I ac io I acustr i ne 
clay  at  the  GB3  inclinometer  location.  The  velocity  gradient  appears  to  be 
fairly  uniform.  Erratic  fluctuations  in  velocity  exist  between  2  and  4  m. 

The  GB3  borehole  description  reveals  that  this  depth  interval  coincides  with 
a  zone  of  small  ice  lenses  that  are  widely  spaced. 

4  Finite  Element  Analyses 

The  stratigraphy  of  the  left  bank  of  the  Great  Bear  River  at  the 
proposed  crossing  site  was  discretized  for  the  finite  element  analysis  as 
shown  in  Figure  4.7.  The  finite  element  grid  consisted  of  22.5  m  of  sand 
overlying  15.5  m  of  medium  to  high  plastic  g I ac io I acustr i ne  clay.  The  stiff 
glacial  till  that  underlies  the  clay  was  not  discretized.  The  finite 
element  grid  was  fixed  at  the  clay/till  interface  to  account  for  the  very 
high  in-situ  modulus  of  the  till. 

The  finite  element  programme  CREEP  computes  movements  in  the  plane  of 
the  slope  only.  Transverse-to-s I  ope  movements  are  assumed  to  be  zero.  The 
inclinometer  data  reveals  that  no  net  transverse-to-s I  ope  movements  were 
recorded  at  hole  locations  GB1A  and  GB3.  At  hole  G82,  the  transverse— to— 
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slope  velocity  recorded  in  the  g I ac io I acustr i ne  clay  is  greater  than  the 
para  I  I e I -to-s I  ope  velocity.  However,  for  reasons  noted  in  the  previous 
section,  the  data  for  hole  GB2  is  not  considered  reliable. 


Elastic  material  properties  were  assigned  on  the  basis 
presented  in  the  literature.  Mean  values  of  bulk  unit  weight  for 
and  clay  were  assigned  on  the  basis  of  measured  values  obtained 
frozen  core  samples.  The  elastic  material  properties  for  the  sand 
are  summarized  in  Table  4.1. 


of  data 
the  sand 
from  the 
and  clay 


TABLE  4.1 

ELASTIC  MATERIAL  PROPERTIES  FOR  FROZEN  GLAC I ODELTA I C 
SAND  AND  GLAC I OLACUSTR I NE  CLAY 


PROPERTY 

SAND 

CLAY 

Young's  Modulus 

Poisson's  Ratio 

Unit  Weight 

4900  MPa 

.495 

1 .84  Mg/m^ 

785  MPa 

.495 

2.04  Mg/ m^ 

Creep  properties  of  the  sand  and  clay  were  altered  to  obtain  close 
agreement  with  the  in-situ  behavior.  Two  finite  element  creep  simulations 
were  carried  out. 
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The  first  creep  analysis  assumed  that  the  sand  and  clay  would  creep  at 
the  same  rate.  A  review  of  the  rheological  behavior  of  frozen  soils 
presented  in  Chapter  I  I  concluded  that  the  stress  versus  strain-rate 
relationship  for  do  I ycrysta II i ne  ice  represents  an  upper  bound  to  the 
constitutive  behavior  of  fine-grained,  ice-rich  permafrost  soils.  Using 
this  argument,  the  flow  law  used  in  the  first  analysis  was: 


r  2.0  x  10~V  (4.1) 

It  was  expected  that  this  flow  law  would  overestimate  the  creep 
velocities  in  the  sand  and  clay.  The  degree  of  overestimation  would  be 
greatest  in  the  sand.  Laboratory  samples  of  the  sand  could  not  be  obtained 
for  laboratory  creep  testing  because  of  sampling  difficulties. 

The  predicted  creep  velocity  profiles  for  all  three  inclinometer 
locations  are  shown  in  Figure  4.8.  The  in-situ  creep  velocities  are  also 
shown  in  this  figure.  At  inclinometer  locations  GB1A  and  GB3,  the  predicted 
velocity  profile  generally  has  a  similar  form  to  the  in-situ  profiles. 
However,  the  magnitude  is  much  greater.  As  expected,  there  is  no 
correlation  between  measured  and  observed  creep  velocities  at  inclinometer 
location  GB2. 

A  close  examination  of  the  magnitude  of  predicted  velocities  versus  the 
observed  values  revealed  that  the  difference  between  the  two  is  a  constant 
factor.  The  creep  strain  rate  varies  linearly  with  the  coefficient  in  the 
simple  power  law.  Hence,  a  proportional  change  in  the  coefficient  will 
reduce  the  predicted  velocities  to  the  observed  values. 
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The  predicted  velocities  were  approximately  six  times  greater  than  the 
in-situ  values.  To  reduce  the  predicted  velocities  proportionate  I y ,  the 
coefficient  in  the  simple  power  law  was  set  equal  to  one  sixth  of  its 
original  value.  Thus,  the  flow  law  for  the  sand  and  clay  now  becomes: 


e  =  .33  x  10  V  (4.2) 
The  exponent  was  left  unchanged  from  its  original  value  of  3.0. 

The  predicted  velocities  reduced  by  a  factor  of  six  are  shown  in  Figure 
4.9.  There  is  good  agreement  between  the  predicted  and  in-situ  velocity 
profiles  for  the  GB1A  and  GB3  inclinometer  locations  through  the  clay 
stratum.  There  is  some  discrepancy  in  the  upper  2  to  3  m  of  the  GB3  hole 
due  to  the  influence  of  the  steel  casing  (Savigny,  1980).  Also,  the  datum 
for  the  finite  element  zero  velocity  at  the  GB1A  and  GB3  holes  were 
translated  horizontally  to  account  for  the  in-situ  movement  which  occurred 
i n  the  under  lying  till. 

The  predicted  in-situ  creep  law  is  shown  in  Figure  4.10  along  with  the 
flow  law  for  ice  (Morgenstern  et  al.,  1980)  and  the  upper  bound  for  soil 
data  as  proposed  by  Nixon  (1978).  The  predicted  flow  law  shows  good 
agreement  with  the  I aboratory  data.  As  shown  in  this  figure,  the  I aboratory 
data  overestimate  the  in-situ  creep  rate. 

The  accumulated  creep  displacement  pattern  after  a  steady  state 
condition  was  achieved  is  shown  in  Figure  4.11.  The  displacement  profiles 
versus  depth  are  indicative  of  non-Newtonian  fluid  behavior. 
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The  preliminary  finite  element  creep  analysis  has  shov/n  fhat  the  i  n- 
situ  creep  behavior  of  the  ice-rich  g I ac io! acustr i ne  clay  can  be  modelled 
with  the  simple  power  law.  However,  the  agreement  between  predicted  and 
observed  velocities  through  the  sand  at  inclinometer  GB1 A  is  less 
convincing.  At  the  outset,  there  was  some  doubt  regarding  the  form  of  the 
constitutive  behavior  of  the  sand  and  whether  the  sand  would  creep  at  all. 
The  ice-rich  g I ac io I acustr i ne  clay  was  expected  to  creep  much  faster  than 
the  overlying  sand.  This  would  set  up  tensile  stresses  in  the  sand  near  the 
sand/clay  interface.  Evidence  of  a  complex  stress  state  at  this  interface 
is  revealed  by  the  observed  erratic  movements  at  the  base  of  the  sand. 

Selection  of  a  suitable  constitutive  relationship  for  the  sand  is 
difficult  due  to  the  lack  of  long  term  low  stress  creep  data  for  natural 
coarse-gra i ned  permafrost  soils.  Creep  tests  on  reconsituted  frozen  sands 
have  shown  that  the  deformations  continue  to  attenuate  with  time  at  stress 
levels  below  the  long  term  strength  (Sayles,  1963).  Creep  deformations  in 
sand  have  also  been  shown  to  be  a  function  of  confining  pressure  (Sayles, 
1973;  Alkire  and  Andersland,  1973).  There  is  also  the  guestion  of  the  form 
of  the  constitutive  relationship  for  the  sand  under  tensile  stresses.  Since 
unfrozen  sand  has  zero  tensile  strength,  all  the  tensile  stresses  must  be 
transmitted  through  the  ice  phase.  Hawkes  and  Me  I  I  or  (1972)  have  shown  that 
the  uniaxial  strength  of  polycrystalline  ice  is  similar  in  tension  and 
compression  at  slow  strain  rates.  However,  the  creep  behavior  of  coarse¬ 
grained  ice-poor  frozen  soils  may  differ  considerably  in  tension  and 
compression.  In  a  compression  creep  test  on  a  coarse-gra i ned  frozen  soil, 
there  are  two  components  which  resist  deformation.  First,  there  are  the 
time-dependent  deformations  occurring  in  the  ice  phase  and  secondly,  inter¬ 
particle  friction  must  also  be  overcome.  In  a  tension  creep  test  only  the 
ice  phase  will  resist  the  application  of  load  since  the  i nterpart i c I e 
friction  has  been  reduced  to  zero.  To  investigate  these  points,  a  second 
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creep  analysis  was  carried  out  in  which  the  sand  was  restricted  to  deform 
elastically.  The  clay  was  assigned  the  flow  law  given  by  equation  4.2. 

The  predicted  velocity  profiles  for  this  analysis  at  all  three 
inclinometer  locations  is  shown  in  Figure  4.12.  Restricting  the  sand  to 
deform  elastically  impedes  the  creep  deformations  in  the  clay  at  the  GB1A 
and  GB2  inclinometer  locations.  At  these  two  inclinometer  locations,  the 
predicted  velocities  do  not  agree  either  in  form  or  magnitude  with  the 
observed  data.  The  predicted  velocities  at  the  GB3  inclinometer  are  only 
slightly  reduced  from  the  analysis  assuming  uniform  creep  properties  because 
this  region  is  well  beyond  the  influence  of  the  sand  stratum. 

The  displacement  pattern  for  the  second  creep  analysis  is  shown  in 
Figure  4.13.  In  this  case,  the  clay  is  being  extruded  out  between  the  sand 
and  underlying  till.  As  the  distance  downslope  from  the  point  the  sand 
pinches  out  increases,  the  displacement  pattern  begins  to  resemble  non- 
Newtonian  fluid  behavior  as  it  did  in  the  first  creep  analysis. 

The  extrusion  of  the  clay  leads  to  very  high  horizontal  tensile 
stresses  in  the  sand  as  shown  in  Figure  4.14.  The  magnitude  of  the  tensile 
stresses  indicate  the  formation  of  tensile  failure  zones  propagating  through 
the  sand  as  it  restricts  the  downslope  movement  of  the  underlying  clay. 

4.5  Assessment  of  Analytical  Results 


On  the  basis  of  the  finite  element  creep  analysis  of  the  left  bank  of 
the  Great  Bear  River  at  the  proposed  Arctic  Gas  crossing,  it  can  be 
concluded  that  the  in-situ  creep  behavior  of  the  g  I  ac io I acustr i ne  clay  can 
be  accurately  represented  by  the  simple  power  law.  The  flow  law  given  by 
equation  4.2  provided  the  best  fif  between  observed  and  predicted 


jjf  I  v  trt 


fc,  9Ut  i '  V  '  '  * 


' 


60 


velocities.  The  exponent  in  the  power  law  is  the  same  as 
polycrystal  line  ice,  however,  the  coefficient  is  one  sixth  of  the 
( Mor gen stern  et  a  I .  1980). 


that  for 
ice  value 


The  in-situ  creep  behavior  of  the  sand  may  also  be  represented  by  the 
simple  power  law.  However,  this  remains  inconclusive  until  further  insight 
is  gained  into  the  in-situ  stress  state  in  the  sand.  Confining  the  sand  to 
behave  elastically  poses  much  too  serious  a  restriction  on  creep  movements 
in  the  clay.  The  in-situ  creep  behavior  indicates  the  clay  creeps  at  a 
faster  rate  than  the  sand.  The  dissimilar  in-situ  creep  behavior  in  the 
clay  and  sand  would  set  up  a  complicated  stress  field  at  the  sand/clay 
interface.  It  is  believed  that  the  faster  creep  movements  in  the  clay  is 
causing  tensile  failure  in  the  sand.  Thus,  to  more  accurately  model  the  f n- 
situ  behavior  of  the  sand,  inclusion  of  a  tensile  failure  criterion  would  be 
necessary.  Approximate  numerical  techniques  are  available  for  treating  soil 
and  rock  as  either  a  'no  tension’  or  ’limited  tension'  material. 
Application  of  such  procedures  is  severely  limited  in  frozen  soil  mechanics 
since  very  little  is  known  about  their  constitutive  behavior  under  tensile 
stress  conditions. 

The  largest  d i screpenc i es  between  observed  and  predicted  velocity 
profiles  appear  to  coincide  with  regions  that  have  experienced  failure.  Two 
such  zones  are  the  base  of  the  sand  at  inclinometer  GB1 A  and  the  entire 
thickness  of  clay  at  GB2  inclinometer  location.  It  is  believed  that  the 
abrupt  change  in  the  in-situ  creep  properties  between  the  clay  and  sand  has 
lead  to  the  development  of  a  zone  of  tensile  failure  in  the  lower  8  m  of 
sand.  The  upper  12  m  of  the  sand,  which  is  well  beyond  the  zone  of 

influence  of  the  sand/clay  interface,  generally  shows  a  uniform  downslope 
velocity. 
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Another  such  zone  occurs  at  the  GB2  inclinometer  location.  The  sand 
deposit  pinches  out  just  upslope  of  this  borehole.  The  stiffer  sand  appears 
to  impede  the  creep  velocity  in  the  underlying  clay.  Downslope  of  the  GB2 
inclinometer,  the  in-situ  creep  behavior  of  the  clay  is  not  influenced  by 
the  overlying  sand  and  creeps  at  a  faster  velocity  than  the  clay  underlying 
the  sand.  The  upslope  restraint  offered  by  the  sand  and  the  downslope  pull 
of  the  unrestrained  creep  in  the  clay  would  lead  to  horizontal  stress 
relaxation  at  GB2.  it  is  believed  that  numerous  slickensided  surfaces 
observed  in  the  core  of  borehole  GB2  occurred  as  a  result  of  continued 
horizontal  stress  relief  while  the  vertical  stress  remained  constant  at  the 
overburden  pressure. 

The  in-situ  creep  behavior  of  the  clay  at  inclinometer  GB3  is  not 
influenced  by  the  sand.  The  steady  state  stress  field  in  the  clay  in  the 
vicinity  of  GB3  lies  well  beyond  the  disturbing  influence  of  the  sand. 

4*6  Modelling  Ice-Rich  Permafrost  as  a  Two-Phase  Continuum 


At  the  proposed  Arctic  Gas  crossing  of  the  Great  Bear  River,  the  data 
gathered  from  inclinometer-  locations  GB1 A  and  GB3  clearly  show  that 
velocities  are  erratic  and  more  movement  is  associated  with  the  large  ice 
lenses  where  they  are  widely  separated  (Savigny,  1980).  The  constant  stress 
creep  tests  carried  out  on  samples  obtained  from  the  field  program  all 
showed  that  the  deformation  behavior  was  affected  by  the  ground  ice 
stucture.  Savigny  (1980)  reported  that  samples  tested  with  a  reticulate  ice 
structure  all  failed  along  the  soil/ice  interface  with  virtually  no  shearing 
through  the  soil.  One  sample  tested  with  stratified  ice  oriented 
perpend icu I ar  to  the  direction  of  load  application  failed  by  radial  tensile 
failure  in  the  soil.  It  is  clear  from  the  foregoing  points,  the  ice 
structures  present  in  an  ice-rich  permafrost  play  an  important  role  in 
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establishing  the  overall  deformation  behavior  of  the  soil  mass.  Thus,  there 
is  good  reason  to  believe  that  the  segregated  ice  features  can  be  treated  as 
discontinuities  in  a  frozen  soil  matrix. 

In  the  field  of  rock  mechanics,  the  presence  of  discontinuities  has 
long  been  recognized  by  engineers  designing  structures  in  the  rock.  The 
behavior  of  a  rock  mass  is  not  only  a  function  of  the  intact  rock  but  also 
the  nature  and  extent  of  the  discontinuities  present. 

Generally,  all  rock  masses  contain  planes  of  potential  weakness  which 
come  in  all  lengths  and  spacings  and  have  varying  degrees  of  influence  on 
the  overal I  rock  mass  properties.  A  geological  field  investigation  can 
accurately  map  the  orientation  and  absolute  position  of  the  more  important 

discontinuities  so  that  their  influence  on  the  rock  mass  behavior  can  be 
adeguately  evaluated. 

convent ional  finite  el emen t  analysis  uti I izes  the  continuum  approach 
which  enforces  material  compatab i I i ty  throughout  the  body  being  analyzed. 
The  continuum  approach  may  be  entirely  adequate  when  dealing  with  small 
joint  spacings  relative  to  the  size  of  the  structure  being  analyzed. 
Situations  where  rocks  deform  along  preexisting  planes  of  weakness  present 
additional  degrees  of  freedom  to  the  rock  mass.  Thus,  the  analytical  model 
must  allow  relative  movement  to  obtain  a  realistic  solution. 

<Vhen  dealing  with  a  single  very  important  discontinuity  whose 
orientation  and  position  are  known  with  confidence,  it  is  possible  to 
represent  plane  of  weakness  explicitly  and  calculate  the  resulting  stresses 
and  deformations  in  the  rock  mass.  Several  interface  elements  have  been 
developed  over  the  past  ten  years  for  solving  problems  that  involve  relative 
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Goodman  et  a  I .  (1968)  introduced  a  four  node  displacement  discontinuity 
element.  Zienkiewicz  et  al.  (1970)  presented  a  continuous  isoparatric 
interface  element  with  linear  strain  along  the  element  and  uniform  strain 
across  its  thickness.  Ghaboussi  et  al.  (1973)  advocates  the  use  of  relative 
displacement  as  an  independent  variable  to  prevent  adjacent  blocks  of 
continuous  elements  from  penetrating  into  each  other.  Simmons  (1981)  used  a 
thin  linear  strain  i soparametr i c  continuum  element  to  model  shearband 
yielding  and  strain  weakening  in  dense,  structured  soils. 

The  essential  feature  in  using  these  interface  elements  is  to  be  able 
to  predict  the  location  of  a  shearzone  or  map  a  through-going  discontinuity. 
Little  is  known  about  the  aerial  extent  of  individual  lenses  of  segregated 
ice  or  reticulate  ice  veins.  The  ice  lenses  may  vary  from  hairline  cracks 
to  massive  ice  bodies  with  vast  horizontal  and  vertical  extent.  The  factors 
governing  the  amount  and  distribution  of  ground  ice  are  the  type  of  host 
material,  availability  of  moisture,  rate  of  freezing  and  geologic  history  of 
events.  The  most  common  form  of  segregated  ice  in  the  g I ac io I acustr i ne  clay 
at  the  proposed  Arctic  Gas  Crossing  site  was  reticulate  ice  with  smaller 
amounts  of  segregated  and  stratified  ice.  The  thicker  primary  veins  in  the 

reticulate  ice  structure  are  typical ly  vertical  making  lateral  correlation 
of  these  ice  features  very  difficult. 

The  variability  of  ground  Ice  is  of  major  significance  when  studying 
the  stability  of  buried,  warm  pipelines  operating  through  permafrost. 
Determination  of  the  amount  and  distribution  of  ground  ice  is  essential  for 
evaluating  the  configuration  of  the  thawed  ditch  bottom  profile.  The 
lateral  variation  of  ground  ice  can  be  determined  by  visual  examination  of 
ditch  walls  or  detailed  logging  of  continuous  undisturbed  samples  obtained 
from  closely  spaced  boreholes. 
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Mapping  of  ground  ice  structures  found  on  ditch  walls  will  provide  the 
most  re  I iable  estimate  of  the  continuity  of  individual  ice  features. 
However,  an  undertaking  of  this  sort  would  be  very  costly.  Ditch  wall 
disturbance  as  a  result  of  overbreak  would  preclude  the  use  of  a  blast  and 
backhoe  operation  to  excavate  the  ditch.  Large  diameter  bucket  wheel 
ditchers  produce  a  smooth  and  relatively  undisturbed  ditch  wall  ideal  for 
observing  the  ground  ice  structure.  Very  few  ditching  trials  have  been 
carried  out  and  reported  in  the  public  literature.  Also,  the  few  trials 
that  have  taken  place  were  primarily  concerned  with  the  ability  of  the 
ditcher  to  excavate  frozen  permafrost  rather  than  mapping  the  ground  ice 
structures . 

Detailed  logging  of  the  ice  structures  observed  in  undisturbed 
permafrost  core  samples  obtained  from  closely  spaced  boreholes  would  be  much 
more  expedient.  The  ground  ice  may  appear  to  be  much  more  variable  in  the 
undisturbed  core  samples  than  that  actually  present  in  the  soil  mass. 
However,  an  estimate  of  the  variability  of  ground  ice  in  the  soil  mass  can 
be  obtained  by  a  statistical  evaluation  of  the  observations  at  individual 
boreho I es . 

Use  of  statistical  analysis  to  solve  variability  problems  in 
geotechnical  engineering  is  relatively  recent.  Quality  control  on 
construction  projects  and  to  a  lesser  extent  the  analysis  of  I aboratory  test 
data  are  two  areas  where  much  of  the  application  of  statistical  techniques 
has  taken  place.  Ideally,  a  known  statistical  model  would  be  used  to 
describe  the  population  under  consideration.  Holtz  and  Kr i zek  (1971) 
studied  several  foundation  case  histories  and  found  that  the  laboratory  test 
data  essentialy  conform  to  the  Gaussian  or  normal  distribution. 
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Speer  et  al.  (1973)  reported  the  results  of  a  ground  ice  variability 
study  undertaken  by  Mackenzie  Valley  Pipeline  Research  Limited  (MVPRL).  The 
purpose  of  the  ice  variability  study  was  to  assess  the  amount  of  ground  ice 
and  its  distribution  in  the  thawed  zone  below  a  warm  oil  pipeline  buried  in 
ice— rich  permafrost.  The  thaw  settlement  potential  of  representative 
geologic  units  through  which  the  propsed  pipeline  might  be  buried  is 
required  to  compute  pipe  stresses  under  operating  conditions. 

When  fine-grained  ice-rich  permafrost  soil  thaws,  water  is  released  and 
settlements  develop  as  the  water  is  expelled  from  the  pore  space  due  to  the 
soils  self  weight  and  applied  I oads .  Laboratory  studies  have  shown  that  for 
light  loadings  in  relatively  ice-rich  soils,  a  major  portion  of  the  thaw 
settlement  occurs  during  the  thaw  stage.  This  component  of  settlement  is 
independent  of  external  loading  and  is  caused  by  excess  water  draining  away. 

A  rough  estimate  of  the  thaw  settlement  for  ice-rich  soils  can  be  obtained 
from  the  thickness  of  visible  ice  lenses.  Thus,  it  would  seem  reasonable  to 
assume  that  the  variation  in  thaw  settlement  provides  an  estimate  of  the 
variability  of  ground  ice. 

Field  and  I aboratory  studies  were  carried  out  by  MVPRL  to  study  ground 
ice  variability  in  silt  and  clay  soils.  The  lateral  variation  of  ground  ice 
content  was  studied  in  rectangular  arrays  of  closely  spaced  borings.  An 
estimate  of  the  total  potential  settlement  was  calculated  from  an  empirical 
correlation  between  frozen  bulk  density  and  thaw  strain  derived  from 
I aboratory  tests  on  undisturbed  core  samples  obtained  during  the  field 
i nvesti gat  ion . 

Two  arrays  were  drilled  in  lacustrine  silts  and  clays  deposited  in  a 
proglacial  lake.  These  deposits  are  considered  textural ly  similar  to  those 
encountered  at  the  proposed  Arctic  Gas  crossing  of  the  Great  Bear  River. 


. 


» 


66 


The  Norman  Wei  Is  study  site  is  located  approximately  3  km  northwest  of  the 
airport  runway  at  Norman  Wei  Is.  The  Landing  Lake  array  is  located 
approximately  53  km  northwest  of  Norman  Wei  Is  on  the  east  bank  of  the 
Mackenzie  River.  At  both  sites,  reticulate  and  segregated  ground  ice 
structures  were  present  in  the  lacustrine  silt  and  clay.  Massive  ice 
structures  were  not  observed. 

The  estimated  total  thaw  settlement  for  the  Norman  Wells  and  Landing 
Lake  study  sites  are  shown  in  Figure  4.15.  The  depth  interval  that  thaw 
settlements  were  calculated  was  2  to  13m.  There  is  good  reason  to  believe 
that  for  a  uniform  stratigraphic  unit,  the  ground  ice  will  be  a  normally 
distributed  variable.  The  combined  data  for  the  Norman  Wells  and  Landing 
Lake  arrays  plotted  in  terms  of  a  thaw  settlement  probability  distribution 
on  normal  probability  paper  is  shown  in  Figure  4.16.  The  observed  data 
shows  little  scatter  about  the  theoretical  normal  distribution.  Also  shown 
in  the  figure  are  the  results  of  the  Ko I mogorov-Sm i rnov  goodness  of  fit  test 
evaluated  for  normality  at  the  20$  significance  level.  This  test  indicates 
that  the  thaw  settlement  in  lacustrine  silts  and  clays  of  the  Mackenzie 
River  Valley  can  be  described  by  a  normal  distribution. 

Although  the  foregoing  has  shown  that  the  variability  of  ground  ice  in 
a  uniform  strat i graph i c  unit  can  be  considered  to  be  normally  distributed, 
it  has  done  I ittle  to  elucidate  the  lateral  extent  of  individual  ice  lenses. 
One  method  of  estimating  the  applicable  horizontal  scale  for  a  given  depth 
interval  is  to  observe  the  surface  profile  of  ground  which  was  once  frozen 
but  has  subsequently  thawed  such  as  seismic  cutlines,  highways  and 
a i rstr i ps . 

A  mathematical  description  of  the  overall  ground  ice  correlation 
between  pairs  of  points  at  various  arbitrary  spacings  is  through  the  concept 
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of  an  auto-correlation  function.  This  function  is  defined  as  the  product  of 
the  deviations  from  the  mean  thaw  settlement  measured  at  pairs  of  points 
spaced  an  arbitrary  distance,  £,  apart  averaged  over  the  entire  population, 
ie: 


n-k 

R(k)  =  J_(  z  (x.-x)(x.  -x)  ] 

n-k  i=i  ' 


(4.3) 


where:  R(k) 
n 

x . 

I 

x 


=  auto-corre I  at  ion  at  spacing  k 
=  total  number  of  observations 

=  observation  at  point. i 
=  mean 


For  the  thaw  settlement  data,  the  auto-correl ation  function  expresses 
the  dependence  of  the  correlation  between  settlements  measured  at  boreholes 
i  and  i  +  k  on  the  distance,  kjj,  between  boreholes.  The  distance  at  which 
the  auto-corre I  at i on  falls  to  zero  indicates  the  horizontal  scale  over  which 
the  thaw  settlements  can  be  correlated.  The  product  of  the  deviations  from 
the  mean  thaw  settlement  spaced  further  than  this  distance  can  be  expected 
to  be  as  often  positive  as  it  is  negative.  Thus,  the  expected  value  of  the 
sum  of  these  products  should  be  zero  beyond  the  correlation  distance  for  a 
very  large  population. 

The  auto-corre I  at  ion  function  determined  from  the  Norman  Wells  and 
Landing  Lake  thaw  settlement  data  is  shown  in  Figure  4.17  The  correlation 
is  almost  zero  for  15  m  and  fluctuates  about  zero  at  greater  spacings 
between  boreholes.  This  suggests  that  thaw  settlements  are  only  weakly 
correlated  at  boreholes  spaced  15  m  apart.  No  data  is  available  at  spacings 
less  than  15  m  since  the  closest  borehole  spacing  at  these  two  ice 
variability  test  sites  was  15  m.  Field  observations  of  disturbed  areas  in 
lacustrine  silts  and  clays  of  the  Mackenzie  River  Valley  indicate  that  15  m 
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represents  an  upperbound  for  the  distance  at  which  ground  ice  content  is 
corre I ated. 

The  thaw  settlement  data  represents  a  one-dimensional  integration  of 
the  ground  ice  conditions  at  a  particular  borehole  over  a  given  depth 
interval.  The  data  derived  from  the  MVPRL  ice  variability  test  sites  at 
Norman  Wells  and  Landing  Lake  suggest  that  for  a  uniform  strat i graph i c  unit, 
the  overall  ground  ice  content  of  the  soil  mass  can  be  described  as  a 
normally  distributed  variable.  The  ground  ice  content  at  a  particular  point 
is  completely  independent  of  neighboring  points  separated  a  distance  of 
15  m.  Thus,  the  time-dependent  deformation  behavior  of  ice-rich  frozen  soil 
is  most  appropr i ate  I y  treated  analytically  by  assuming  homogeneity  for  the 
ice-rich  permafrost  until  a  better  understanding  of  the  lateral  extent  of 
individual  ice  lenses  is  obtained. 
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FIGURE  4.1  SITE  PLAN  OF  PROPOSED  ARCTIC  GAS  CROSSING  OF 

GREY  BEAR  RIVER,  N.W.T. 

(Savigny,  1980) 
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FIGURE  4.2 


SITE  PLAN  OF  LEFT  BANK  OF  GREAT  BEAR  RIVER  AT 
PROPOSED  ARCTIC  GAS  CROSSING 
(Savigny,  1980) 
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FIGURE  4.3  STRATIGRAPHIC  CROSS  SECTION;  LEFT  BANK  OF  GREAT  BEAR 

RIVER  AT  PROPOSED  ARCTIC  GAS  CROSSING 
(Savigny,  1980) 
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FIGURE  4.4  THERMAL  CROSS  SECTION;  LEFT  BANK  OF  GREAT  BEAR 

RIVER  AT  PROPOSED  ARCTIC  GAS  CROSSING 
(Savigny,  1980) 
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FIGURE  4.5 


SITE  PLAN  OF  INCLINOMETER  LOCATIONS  ON  LEFT  BANK  OF 
GREAT  BEAR  RIVER  AT  PROPOSED  ARCTIC  GAS  CROSSING 
(Savigny,  1980) 
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FIGURE  4.6  VELOCITY  VERSUS  DEPTH  PROFILES  FOR  GBIA,  GB2,  AND  GB3 
INCLINOMETER  LOCATIONS 
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GREAT  BEAR  RIVER 


FIGURE  4.7  FINITE  ELEMENT  DISCRETIZATION  OF  PROPOSED  ARCTIC 

GAS  CROSSING  OF  THE  GREAT  BEAR  RIVER 
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FIGURE  4.8 


PREDICTED  STEADY  STATE  CREEP  VELOCITY  PROFILES, 
FIRST  CREEP  ANALYSIS 
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FIGURE  4.9  PREDICTED  STEADY  STATE  VELOCITY  PROFILES,  FACTORED 

VELOCITIES  FROM  FIRST  CREEP  ANALYSIS 


PREDICTED  FLOW  LAW  FOR  GLACIOLACUSTR INE  CLAY 
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FIGURE  4.10  COMPARISON  OF  IN-SITU  CREEP  LAW  AND  STEADY-STATE  CREEP  DATA  FOR 

ICE-RICH  GLACIOLACUSTR  INE  CLAY  FROM  VARIOUS  LOCATIONS  IN  THE 
MACKENZIE  RIVER  VALLEY,  N.W.T. 


Deformations  Exaggerated  250  x 
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FIGURE  4.11  DEFORMED  GRID  FOR  STEADY  STATE  CREEP  IN 

SAND  AND  CLAY 
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FIGURE  4.12  PREDICTED  STEADY  STATE  VELOCITY  PROFILES, 

CREEP  IN  GLACIOLACUSTRINE  CLAY  ONLY 
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FIGURE  4.13  DEFORMED  GRID  FOR  STEADY  STATE  CREEP  IN  CLAY  ONLY 
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FIGURE  4.14  HORIZONTAL  STRESS  CONTOURS  FOR  STEADY  STATE 

CREEP  IN  CLAY  ONLY 
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THAW  SETTLEMENT  DATA  FOR  NORMAN  WELLS  AND  LANDING  LAKE 

ICE  VARIABILITY  TEST  SITES  (Speer  et.  al„  1973) 
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THEORETICAL  NORMAL 
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FIGURE  4.16  PROBABILITY  RELATIONSHIP  FOR  THAW  SETTLEMENT 

IN  ICE-RICH  LACUSTRINE  SILT  AND  CLAY  (Speer  et.  al.,  1973) 
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AUTO  -  CORRELATION  FUNCTION  FOR  THAW 
SETTLEMENT  IN  ICE  -  RICH  LACUSTRINE  SILT 
AND  CLAY 


CHAPTER  V 


IN-SITU  CREEP  ANALYSIS  OF  FOX  TUNNEL 


5.1  General 


The  Fox  Tunnel  is  located  approximately  15  krn  north  of  Fairbanks, 
Alaska.  The  tunnel  was  excavated  into  perenially  frozen  silts.  These  silts 
have  been  of  an  economic  interest  in  the  Fairbanks  area  since  the  turn  of 
the  century  because  of  the  need  to  remove  thick  sections  of  the  material  to 
expose  the  underlying  gold-bearing  gravels.  The  old  mining  practice 
consisted  of  stripping  the  overburden  of  perenially  frozen  silt  to  expose 
the  gravels  which  were  subsequently  worked  with  dredges.  The  Fox  Tunnel 
excavations  were  part  of  two  separate  research  programs  from  1966  to  1969. 
The  objectives  of  these  research  projects  were  to  investigate  methods  of 
subsurface  exploration  in  perenially  frozen  ground  and  evaluate  subsurface 
openings  as  a  shelter,  storage  space  and  site  for  military  activities. 

The  first  stage  of  this  project  was  to  excavate  a  tunnel  110  m  long  in 
the  winters  from  1963  to  1966  as  shown  in  Figure  5.1.  The  excavation  was 
carried  out  using  a  continuous  mining  machine  and  a  modified  blasting 
technique.  The  tunnel  was  excavated  through  ice-rich  Fairbanks  silt  for  its 
entire  length.  The  tunnel  portal  was  excavated  into  a  near-vert ica I  silt 
escarpment  left  from  old  placer  mining  operatidns.  Closure  measurements 
were  undertaken  following  the  mining  operations.  In  November,  1966,  a 
1.22  m  ventilation  shaft  was  drilled  from  the  surface  and  connected  to  the 
end  of  the  tunnel.  The  ventilation  shaft  was  left  open  the  first  winter  to 
allow  the  heavier  cold  winter  air  to  circulate  through  the  tunnel  by  natural 
convection.  The  cooler  air  lowered  the  tunnel  wall  temperature  from  -1.0°C 
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to  -5.7°C  (McAnerney,  1968).  The  lower  ground  temperature  substantially 
decreased  the  closure  rate  of  the  tunnel.  Deformations  have  been  negligible 
since  the  ventilation  stack  was  drilled. 

The  second  stage  of  this  project  involved  excavating  a  winze  from  the 
existing  tunnel  to  provide  access  to  the  strati  graph ica I  I y  lower  gold 
bearing  gravels  and  bedrock.  This  phase  of  the  project  is  shown  as  the 
shaded  portions  in  Figure  5.2.  This  phase  was  a  cooperative  effort  between 
the  U.S.  Bureau  of  Mines  (USBM)  and  the  U.S.  Army  Cold  Regions  Research  and 
Engineering  Laboratory  (USA  CRREL) .  The  winze  terminates  in  a  room; 
designated  as  Room  B  by  the  USA  CRREL.  A  second  room  USBM  Room  A,  was 
excavated  south  of  the  inclined  winze.  The  winze  and  two  rooms  were 
excavated  using  air  and  conventional  blasting  techniques.  Slabbing  of  the 
gravel  occurred  in  regions  where  the  thickness  of  gravel  separating  the 
ceiling  from  the  silt  was  thin.  This  necessitated  that  this  material  be 
scaled  from  the  roof  to  ensure  a  safe  working  area.  Extensive  deformation 
studies  were  carried  out  in  both  rooms  (Thompson  and  Sayles,  1972; 
Pettibone,  1973).  These  will  be  discussed  in  a  subsequent  section. 

5.2  Site  Geology 

The  sequence  of  geologic  sediments  encountered  in  the  main  tunnel  and 
winze  excavation  exposures  are  shown  in  Figure  5.3.  The  disintegrated 
schist  bedrock  at  the  tunnel  site  is  mantled  by  unconsolidated,  ice-cemented 
silts  and  gravels.  Immediately  overlying  the  bedrock  are  the  gold-bearing 
gravel  deposits.  The  gravel  is  capped  by  a  thick  re- transported  silt 
depos i t. 

The  floor  of  both  rooms  excavated  at  the  end  of  the  winze  consist  of 
Birch  Creek  schist.  The  bedrock  is  perenially  frozen  and  deeply  weathered 
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in  ■‘"hese  exposures.  X-ray  diffraction  analyses  performed  on  samples  of  the 
decomposed  bedrock  just  below  the  grave  I /bedrock  interface  showed  a  high 
montmor i I  I  on i te  peak  (Sellman,  1972). 

A  stream  deposited  gravel  layer  lies  unconformab I y  on  top  of  the 
bedrock.  Stream  deposition  was  inferred  from  the  strong  imbrication  of  the 
pebbles  and  cobbles.  The  average  thickness  of  the  gravel  deposit  in  this 
region  is  approximately  4  m.  Pockets  of  fine-grained  materials  are 
interbedded  in  the  sand.  The  gravel  is  perenially  frozen  and  bonded  with 
ice.  Pore  ice  is  visible  in  the  voids  but  the  particles  retain  grain-to- 
gra i n  contact. 

A  thick  deposit  of  silt,  the  uppermost  unit  at  this  site,  lies 
unconformab I y  on  top  of  the  gravel.  Although  irregular,  the  silt-gravel 
contact  is  sharp.  The  silt  deposit  ranges  in  age  from  I  I  I i no i an  to  recent. 
The  lllinoian  deposits  contain  considerably  less  ground  ice  and  organic 
matter  than  the  younger  units.  The  overlying  silt  of  Wisconsin  age  is 
character i zed  by  large  ice  wedges.  Retransported  silts  of  recent  age  mantle 
the  top  of  the  upland  plain  at  the  tunnel  site.  Generally,  there  are  two 
types  of  ground  ice  present  in  the  frozen  silt;  ice  lenses  and  masses  formed 
by  ice  segregation,  and  the  massive  ice  structures  of  Aufeis  (buried  ice). 
Ice  volumes  range  from  54  to  79$. 

5.3  In-Situ  Deformation  Studies 


The  in-situ  creep  movements  of  the  ice-rich  Fairbanks  silt  were  studied 
independently  by  Swinzow  (1970),  Thompson  and  Sayles  (1972)  and  Pettibone 
(1973) . 
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The  closure  rate  of  the  main  tunnel  was  measured  by  Swinzow  (1970) 
immediately  following  the  final  excavation  season.  No  details  on  the  type 
of  instrumentation  or  the  exact  time  the  measurements  were  recorded  are 
given.  Horizontal  and  vertical  closure  measurements  were  recorded  at  two 
different  locations,  Sta.  1+50  and  Sta.  3+50  (45.7  m  and  106.7  m, 
respectively,  from  the  portal).  The  closure  measurements  recorded  at  the 
two  tunnel  sections  are  shown  in  Figures  5.4.  The  air  temperature  in  the 
tunnel  at  the  time  the  readings  were  recorded  was  very  close  to  the  melting 
point  of  ice  (McAnerney,  1968).  Measurements  recorded  over  a  3200  hour  time 
span  showed  that  the  deformation  rate  continued  to  attenuate  with  time. 
Very  little  difference  was  noted  between  horizontal  and  vertical  closure  of 
the  tunnel  at  both  locations.  A  closure  of  14.6  cm  was  recorded  at  the 
instrumented  section  furthest  removed  from  the  tunnel  portal.  The 
difference  in  overburden  between  these  two  sections  was  only  2.1  m.  The 
difference  in  closure  between  the  two  measuring  stations  is  likely  due  to 
temperature  since  the  end  of  the  tunnel  was  consistently  warmer  during  the 
observation  period  (Swinzow,  1970).  The  closure  measurements  clearly  showed 
that  artificial  ventilation  was  inadequate  to  control  the  deformations 
within  tolerable  limits.  In  December,  1965  a  vertical  ventilation  shaft  was 
drilled  at  the  end  of  the  tunnel.  Cold  winter  air  created  by  the  chimney 
effect  lowered  the  temperature  of  the  surrounding  silt  to  -6.7°C  and  -10°C 
at  distances  of  0.61  m  and  2.35  m  from  the  tunnel  wall.  The  colder  tunnel 
wall  temperature  has  subsequently  reduced  closure  substantially. 

McAnerney  (1968)  measured  the  vertical  deformation  of  the  tunnel  at 
eight  deformation  points  as  shown  in  Figure  5.5(a).  The  readings  plotted 
for  one  year  commencing  in  April,  1966,  are  shown  in  Figure  5.5(b).  The 
readings  recorded  by  Swinzow  in  1965  are  superimposed  on  this  plot.  Note 
that  the  deformation  is  negligible  after  the  ventilation  shaft  is  opened  and 
cold  winter  air  enters  the  tunnel  by  natural  convection. 
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Field  observations  of  the  USBM  room  were  reported  by  Pettibone  (1973). 
The  USBM  room  was  initially  excavated  4.6  x  15.2  x  2.4  m  high.  In 
subsequent  weeks  2.3  m  was  slabbed  from  each  side  to  widen  the  room  to 
approximately  9  m.  The  room  was  instrumented  with  vertical  deformation 
points  as  excavation  proceeded. 

The  instrumentation  was  designed  to  measure  mass  movement  independent 
from  parting.  Floating  points  used  to  measure  mass  subsidence  consisted  of 
1.8  m  long  steel  rock  bolts  embedded  in  the  roof  in  overbored  holes.  the 
parting  was  measured  with  a  1 52  mm  wood  dowel  placed  adjacent  to  the  rock 
bolt.  A  0.6  mm  long  rock  bolt  was  permacreted  into  the  bedrock  floor 
directly  underneath  the  rock  bolt  and  wood  dowel.  Two  thermocouple  strings 
were  installed  in  the  roof  and  one  in  the  wall. 

The  in-situ  creep  deformation  versus  time  plot  for  one  of  the  measuring 
points  is  shown  in  Figure  5.6.  Deformation  readings  commenced  in  early 
February,  1969.  Subsidence  of  the  USBM  room  increased  rapidly  as  the  room 
was  enlarged.  Following  excavation,  the  deformation  rate  dropped 
substantially.  The  room  temperature  up  to  mid-April  was  approximately 
-6.7°C.  From  mid-April  to  mid-July,  1969,  the  room  temperature  increased  to 
about  -3.3°C.  The  accompanying  increase  in  deformation  rate  caused  some 
thin  slabbing  near  the  back  of  the  room.  In  early  August,  1969,  an 
insulated  bulkhead  was  installed  at  Sta.  1t20.  Following  this  the 
temperature  increased  rapidly  and  was  maintained  at  approximately  -2.2°C. 
The  bulkhead  was  removed  in  early  January,  1970.  During  the  period  the 
bulkhead  was  in  place,  the  rate  of  floor  to  roof  closure  continually 
exceeded  the  closure  between  the  floor  and  1.8  m  depth  into  the  roof. 
Hydraulic  props,  installed  as  a  safety  precaution  under  portions  of  the  roof 
that  were  scaling,  continued  to  pick  up  load  during  this  period  indicating 
that  they  were  picking  up  load  from  overall  mass  subsidence  in  addition  to 
the  dead  load  of  the  slab. 


,  -  i  -  *  !  7  :;»  r.fifi  V  ■  '•£’*  “  '  ' 


91 


After  the  bulkhead  was  removed,  the  closure  rate  has  decreased  to 
approximately  1.5  cm/yr.  This  is  consistent  with  the  observations  in  the 
CRREL  room.  Colder  average  temperatures  for  the  duration  of  observations  in 
the  USBM  room  resulted  in  less  closure  in  this  room  than  the  CRREL  room.  In 
all,  30  cm  of  floor  to  roof  deformation  was  recorded  in  37  months  with 
parting  accounting  for  approximately  20$  of  total  room  closure.  Several 
cracks  were  observed  in  the  roof  in  April,  1970.  Pettibone  (1973)  concludes 
that  vertical  separation  and  expansion  of  the  formation  account  for  the 
measured  parting. 

Thompson  and  Sayles  (1972)  reported  the  results  of  a  finite  element 
simulation  of  the  closure  of  the  USA  CRREL  room  in  the  Fox  Tunnel. 
Horizontal  and  vertical  closure  measurements  were  reported  for  one  year. 
Laboratory  samples  were  obtained  from  uniaxial  creep  testing. 

The  authors  used  the  simple  power  law  to  predict  the  in-situ  behavior 
of  the  room.  They  concluded  that  the  form  of  the  flow  law  predicted  the  in- 
situ  behavior  guite  well.  However,  the  back-calculated  flow  law  predicted  a 
creep  rate  3.3  times  faster  than  had  been  established  in  the  laboratory  on 
the  same  frozen  silt.  A  few  comments  can  be  made  regarding  the  analysis 
carried  out  by  Thompson  and  Sayles  (1972). 

A  two  dimensional  finite  element  analysis  was  used  to  model  the  room. 
The  actual  field  configuration  was  three-dimensional.  It  is  common  practice 
to  assume  that  sections  within  three  tunnel  diameters  will  be  influenced  by 
the  end-restra i nt  offered  by  the  tunnel  face.  Both  of  the  instrumented 
sections  in  this  room  are  within  this  distance.  The  end  effects  of  the  room 
would  delay  the  onset  of  in-situ  steady-state  creep. 
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There  is  a  slight  error  involved  in  the  boundary  condition  assumed 
along  the  grave l-bedrock  interface.  In  their  analysis,  Thompson  and  Sayles 
(1972)  assume  the  base  of  the  gravel  is  free  to  translate  in  the  horizontal 
direction  due  to  a  lean  clay  that  is  present  immediately  below  the  gravel- 
bedrock  contact.  The  decomposed  bedrock  exposed  at  the  base  of  the  gravels 
in  this  room  are  perenially  frozen  and  no  excess  segregated  ice  features 
were  observed  (Sellman,  1972).  Thus,  the  clay  would  offer  some  frictional 
resistance  which  would  decrease  the  predicted  creep  rate. 

The  form  of  the  constitutive  relationship  of  the  placer  gravels  is 
questionable.  The  initial  finite  element  simulations  indicated  that  in-situ 
behavior  of  the  room  could  not  be  predicted  assuming  constitutive  models 
which  exhibit  continuously  decreasing  creep  rate  with  time  (Thompson  and 
Sayles,  1972).  The  simple  power  law  for  the  gravel  was  arrived  at  by  a  best 
fit  to  the  actual  field  curves.  The  average  density  and  moisture  content  of 
the  gravel  is  2.08  Mg/m3  ancj  10%,  respectively.  Sellman  (1972)  described 
the  gravel  as  retaining  gra i n-to-gra i n  contact  with  no  massive  ice  forms 
present. 

In  light  of  the  review  of  the  constitutive  behavior  of  frozen  soils 
presented  in  Chapter  2,  it  would  appear  that  the  behavior  of  the  frozen 
gravel  will  be  governed  by  frictional  response.  Swinzow  (1970)  reported 
that  a  room  excavated  in  cold,  densly  packed,  permafrost  till  did  not  deform 
in  excess  of  measurement  error  after  3  years  of  operation.  Thus,  depending 
on  temperature,  the  in-situ  deformation  behavior  of  the  gravel  may  be  more 
accurately  represented  by  elastic  behavior  or  a  creep  law  which  is  a 
function  of  confining  pressure. 
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Examination  of  the  laboratory  creep  data  of  the  Fairbanks  silt  is 
necessary  to  establish  whether  a  steady  state  condition  did  exist. 
Unconfined  compression  creep  tests  were  performed  on  undisturbed  samples  of 
frozen  Fairbanks  silt  obtained  from  the  tunnel  exposure.  The  test 
temperature  was  -1.7°C  and  the  stress  levels  ranged  from  245  kPa  to  2100 
kPa.  A  simple  power  law  was  fitted  to  the  I  aboratory  creep  data  with  the 
following  parameters: 


e  =  1.16  x  10  \  4  (5.1) 
where  the  units  for  strain  rate  and  stress  are  yeai  •  and  kPa,  respectively. 

Thompson  and  Sayles  (1972)  reported  strains  in  the  I aboratory  typically 
in  excess  of  40#  in  8  days.  Tests  performed  by  Savigny  (1980)  on  ice-rich 
fine-grained  g I ac io I acustr i ne  clay  showed  that  compatab i I i ty  at  the  soil-ice 
interface  would  not  exist  at  strains  of  this  magnitude.  Thus,  the  minimum 
strain  rates  measured  in  the  laboratory  may  well  exceed  the  steady-state 
rate  and  the  exponent  of  4  is  considered  high. 

A  more  detailed  investigation  of  the  U.S.  Army  CRREL  Room  was 
undertaken  to  gain  further  insight  into  the  in-situ  creep  behavior  of  the 
frozen  silt. 

The  horizontal  and  vertical  closure  of  the  USA  CRREL  room  were 
monitored  at  two  instrumented  test  sections,  Sta.  1+83  and  Sta.  2+02.  The 
horizontal  and  vertical  room  deformations  at  each  section  were  recorded 
between  four  closure  points  installed  in  the  roof,  walls  and  floor.  Each 
point  consisted  of  a  32  mm  diameter  wooden  dowel  embedded  0.3  m  into  the 
permafrost.  A  steel  tape  in  conjunction  with  a  steel  ruler  with  an  accuracy 
of  1.6  mm  (1/16  inch)  was  used  to  measure  the  closure  between  opposing 
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points  in  the  walls,  and  roof  and  floor.  In-situ  ground  temperatures  were 
monitored  by  four  strings  of  thermocouples  installed  at  Sta.  1+96.  The 
thermocouples  were  installed  in  the  roof,  floor  and  each  side  wall.  The 
ground  temperatures  were  recorded  at  0.3  m  intervals  from  the  walls  and  roof 
to  a  total  depth  of  1.83  m. 

The  closure  and  temperature  instrumentation  could  not  be  installed  for 
three  weeks  after  the  room  was  excavated.  An  insulated  bulkhead  was 
installed  at  the  entrance  to  the  room  (Sta.  1+81)  approximately  eight  weeks 
after  the  initial  closure  and  temperature  readings  were  recorded.  The 
purpose  of  this  bulkhead  was  to  maintain  a  constant  room  temperature  thus 
minimizing  the  effect  of  seasonal  fluctuations  in  air  temperature  on  the 
permafrost  temperature.  Approximately  18  weeks  later,  a  second  bulkhead  was 
installed  by  the  USBM  at  Sta.  1+20  and  subsequently  removed  five  months 
after  i nsta I  I  at  ion . 

The  vertical  and  horizontal  closure  measurement  recorded  at  Sta.  1+83 
and  2+02  in  the  USA  CRREL  room  are  shown  in  Figures  5.7  and  5.8, 
respectively.  The  in-situ  ground  temperatures  measured  in  the  roof  and 
walls  are  also  plotted  in  Figures  5.7  and  5.8.  The  vertical  and  horizontal 
closure  rate  determined  from  the  closure  versus  time  data  (Figures  5.7  and 
5.8)  is  shown  in  Figure  5.9.  The  closure  observations  were  recorded  for  a 
period  of  3.65  years. 

The  vertical  and  horizontal  closure  measurements  increased  rapidly 
during  the  first  year  of  observations.  This  is  probably  due  to  the  general 
warming  trend  in  the  surrounding  permfrost  as  indicated  by  the  temperature 
measurements.  After  removal  of  the  USBM  bulkhead  at  Sta.  1+20,  the  rate  of 
closure  decreased  as  a  consequence  of  the  cooler  permafrost  temperatures . 
The  greatest  vertical  and  horizontal  movements  were  recorded  at  Sta.  2+02. 


' 


% 


‘ 


„ 


95 


After  3  1/2  years  of  observations,  48  cm  and  16.5  cm  of  vertical  and 
horizontal  movement,  respectively,  had  taken  place. 

The  temperature  dependence  on  in-situ  creep  behavior  is  clearly 
illustrated  by  the  vertical  closure  at  Sta.  1+83.  The  first  bulkhead  was 
installed  within  .75  m  of  this  instrumented  section.  The  colder  ambient  air 
temperature  just  outside  of  the  room  would  lower  the  permafrost  temperature 
at  Sta.  1+83.  This  would  explain  the  slower  vertical  deformation  rate  at 
Sta.  1+83.  From  January,  1970  to  October,  1972,  the  deformation  readings 
indicate  a  steady-state  closure  rate  at  Sta.  1+83  at  a  velocity  of  1  cm/yr. 
For  the  duration  of  this  time  interval  the  permafrost  temperature  has  been 
maintained  at  an  average  value  of  -2.2°C. 

The  vertical  closure  at  Sta.  2+02  shows  much  less  sensitivity  to 
changes  in  permafrost  temperature  brought  about  by  installation  and  removal 
of  the  USBM  bulkhead.  At  this  section,  the  closure  rate  continued  to 

attenuate  with  time  for  the  entire  period  of  time  closure  obsevations  were 
recorded.  During  the  final  ten  months  of  observations,  the  closure  rate  at 
Sta.  2+02  was  2.5  cm/yr. 

The  horizontal  closure  at  Sta.  1+83  and  2+02  continued  to  accelerate 
for  the  first  year  movements  were  recorded.  Following  removal  of  the  USBM 
bulkhead,  the  horizontal  movements  at  Sta.  1+83  fluctuated  seasonally  with 
1.3  cm  of  net  closure  being  recorded  for  the  final  2.5  years  of 
observations.  During  the  warm  summer  months  and  early  fall,  1  to  2  cm  of 
horizontal  movement  was  recorded.  Very  little  or  negative  closure  was 
recorded  during  the  cold  winter  months.  Attenuation  of  the  horizontal 
closure  rate  at  Sta.  2+02  was  observed  for  1  year  following  removal  of  the 
USBM  bulkhead.  Seasonal  fluctations  have  subsequently  taken  place. 
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centimetre  of  net  horizontal  closure  was  recorded  during  the  final  1.7  years 
of  obsevations. 

A  finite  element  simulation  of  the  USA  CRREL  room  was  carried  out  using 
the  programme  CREEP.  The  purpose  of  this  simulation  was  to  assess  the 
validity  of  the  simple  power  low  for  ice  as  an  upper  bound  solution  for  the 
in-situ  creep  deformations  recorded  in  the  frozen  silt  while  the  frozen 
p I acer  gravels  were  treated  as  an  elastic  material. 

The  discretized  continuum  of  the  CRREL  room  used  for  this  finite 
element  simulation  is  shown  in  Figure  5.10.  Only  the  Fairbanks  silt  and 
placer  gravels  were  discretized.  It  was  assumed  that  the  competent  bedrock 
floor  would  not  contribute  to  vertical  closure  of  the  room.  Only  one  half 
of  the  tunnel  was  analyzed  due  to  symmetry.  The  thickness  of  overburden 
above  the  roof  center  I ine  was  16.9  m.  The  ground  surface  was  assumed  stress 
free.  The  two  vertical  boundaries  were  fixed  by  roller  constraints  which 

allowed  movement  in  the  vertical  direction  only.  The  base  of  the  gravel  was 
placed  on  roller  constraints  which  allowed  movement  in  the  horizontal 
direction  only  as  was  done  in  finite  element  analysis  carried  out  by 
Thompson  and  Sayles  (1972).  The  elastic  material  properties  of  the  frozen 
silt  and  gravel  are  given  in  Table  5.1. 

The  Young's  Modulus  values  and  unit  weights  were  obtained  from  Thompson 
(1970).  Creep  properties  of  the  silt  and  gravel  were  assigned  the  values  of 
pure  ice  at  -2°C  (Morgenstern  et  al.,  1980),  i.e.: 

_8  3 

e  =  2.0  x  10  a  (5.2) 


where  the  units  of  stress  and  strain  rate  are  kPa  and  yr  ^ ,  respectively. 
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TABLE  5.1 

ELASTIC  MATERIAL  PROPERTIES  FOR 
FROZEN  SILT  AND  GRAVEL 
FOX  TUNNEL,  ALASKA 


PROPERTY 

SILT 

GRAVEL 

Young's  Modulus 

346  MPa 

713  MPa 

Poissons'  Ratio 

.495 

.450 

Unit  Weight 

1  .47  Mg/m3 

2.08  Mg/m3 

The  nodal  point  label  led  V  in  Figure  5.10  was  chosen  to  represent  the 
in-situ  vertical  closure  points.  The  vertical  displacement  of  node  V  was 
assumed  to  represent  the  total  vertical  closure  since  the  closure  point 
installed  in  the  floor  was  fixed  into  bedrock. 

The  vertical  steady-state  velocity  predicted  by  the  programme  CREEP 
using  the  simple  power  law  given  by  equation  5.2  to  model  the  creep  behavior 
of  the  frozen  silt  was  1  mm/yr.  The  accummulated  closure  of  the  node  point 
labelled  H  in  Figure  5.10  after  2.5  years  was  11  mm.  The  horizontal 
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displacement  of  node  H  was  assumed  to  be  one  half  of  the  total  hoizontal 
closure.  The  predicted  results  are  at  least  one  order  of  magnitude  less 
than  the  in-situ  measurements.  The  design  of  the  instrumentation  installed 
in  the  CRREL  room  was  such  that  vertical  separation  in  the  ceiling  could  not 
be  measured  independently  of  total  floor  to  ceiling  closure.  In  the  USBM 
room,  20 %  of  the  total  room  closure  was  due  to  parting  in  the  ceiling. 
Thus,  a  portion  of  the  total  closure  measured  in  the  CRREL  room  is  likely 
due  to  vertical  separation  in  the  ceiling. 

In  this  case,  the  simple  power  law  for  ice  does  not  represent  a  valid 
upper  bound  for  the  in-situ  creep  behavior  of  the  frozen  silt.  However, 
examination  of  the  closure  rate  versus  time  data  presented  in  Figure  5.9 
clearly  shows  that  steady-state  conditions  did  not  exist  during  the  first 
year  of  observations  as  assumed  by  Thompson  and  Sayles  (1972).  It  is 
entirely  possible  that  the  overall  closure  measurements  that  were  recorded 
consisted  of  both  plastic  and  creep  flow.  Plastic  flow  in  this  sense  refers 
to  time-dependent  failure  of  the  frozen  soil. 

5*3  Assessment  of  In-Situ  Deformation  Behavior 


The  in-situ  deformation  studies  of  the  Fox  Tunnel  reported  by  Swinzow 
(1970),  Thompson  and  Sayles  (1972)  and  Pettibone  (1973)  all  showed  that  the 
deformations  in  the  ice-rich  Fairbanks  silt  became  excessive  when  the  room 
air  temperature  approached  the  ambient  soil  temperature  of  approximately 
-1.7°C.  In  each  case,  the  room  temperature  was  lowered  to  restrain  the 
movements  within  tolerable  limits. 

The  stress  release  in  the  ceiling  of  the  USA  CRREL  room  and  the  USBM 
room  initiated  by  excavating  the  underlying  frozen  placer  gravels  is  similar 
to  that  occurring  in  a  soil  mass  located  above  a  yielding  trap  door.  In  the 


' 


'"•x  '  '  -  ‘  -i  •  ?  W 


> 


* 


99 


past,  it  has  been  recognized  that  the  behavior  of  a  soil  immediately  above  a 
flexible  underground  structure  and  that  of  a  granular  mass  above  a  yielding 
bottom  of  a  silo  are  equivalent  for  all  practical  purposes  (Ladanyi  and 
Hoyaux,  1969).  This  assumption  allows  arching  or  bin  theory  to  be  used  to 
estimate  the  contact  pressures  between  the  ground  and  an  underground 
structure. 

The  general  expression  for  the  pressure  acting  on  the  roof  of  a  tunnel 
derived  from  arching  theory  is  given  by: 

a,  =  B(y-c/B)  ( 1  _e""K+an4>D/B)  (5.3) 

Ktan<j> 

where:  =  vertical  contact  pressure 

B  =  half  width  of  the  zone  of  arching 

(half  width  of  the  underground  rooms  in  this  case) 

D  =  depth  of  the  roof  below  the  surface 
Y  =  unit  weight  of  soil  above  the  roof 
c  +  4)  =  shear  strength  parameters  of  the  soi  I  above  the  roof 
K  =  ratio  between  vertical  and  horizontal  stress 

It  is  clear  from  the  above  expression  that  an  underground  opening 
located  in  a  cohesionless  material  will  not  remain  stable  without  additional 
ground  support  provided  by  a  tunnel  liner  or  timber  lagging.  For  an  unlined 
opening  in  a  cohesive  material,  the  following  relationship  must  be 
sat i sf i ed : 


C  >  By 


(5.4) 
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If  the  cohesion  is  less  than  Bt ,  the  roof  of  the  underground  opening 
must  be  supported.  For  the  USA  CRREL  room,  values  of  3  =  2.44  m  and  Y  = 

1.47  Mg/m-^  require  a  minimum  cohesion  of  35.1  kPa  to  ensure  stability  of  the 
room. 


Vialov  (1962)  proposed  that  the  strength  of  a  frozen  soil  can  be 
represented  by  a  series  of  Mohr-Coulomb  failure  envelopes,  each  one 
correspond i ng  to  a  specific  time  to  failure.  The  relationship  between 
maximum  shear  and  normal  stress  takes  the  form: 

t(+)  =  c+  +  an+an<t,+  (5.5) 


where:  x(t) 


a, 


v  *+ 


shear  strength  at  time  =  t 
maximum  normal  stress  on  shear  plane 
cohesion  and  angle  of  internal  friction,  respectively 
(functions  of  time  and  temperature) 


The  time  variation  of  the  cohesion  intercept,  c^.,  can  be  described  by 

c f  =  6/ I n ( t  f/B)  (5.6) 


where:  B,  B  =  experimentally  determined  material  parameters 

=  t ime  to  f a i I ure 

Time  and  temperature  have  very  little  influence  on  the  angle  of  internal 
friction. 

The  laboratory  studies  carried  out  to  date  indicate  that  the  long  term 
strength  of  frozen  soil  is  frictional.  The  angle  of  internal  friction  can 
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well  be  approximated  by  the  effective  friction  angle  for  the  same  soil  in  an 
unfrozen  state.  For  coarse-gra i ned  soils,  the  friction  angle  is  relatively 
independent  of  temperature  and  strain  rate.  Frictional  resistance  is 
hindered  by  the  presence  of  ice  in  ice-rich  frozen  soils.  In  very  ice-rich 
soils  which  have  very  little  gra i n-to-gra i n  contact,  the  strength  will  be 
that  of  the  ice  and  in  the  long  term  will  approach  zero. 

The  general  nature  of  the  failure  envelopes  for  frozen  soils  is 
nonlinear.  The  time  variation  of  the  cohesion  and  angle  of  internal 
friction  is  reflected  by  the  mutual  distribution  of  each  envelope.  The 
inclination  of  each  envelope  represents  the  time  variation  of  the  angle  of 
internal  friction.  The  loss  of  cohesion  with  time  is  represented  by  the 
cohesion  intercept  of  each  envelope. 

In  predicting  the  time-dependent  settlements  for  foundations  in  frozen 
soils,  Ladanyi  (1981)  recognized  that  creep  and  consolidation  can  occur 
simultaneously  in  frozen  soils  containing  large  amounts  of  unfrozen  water. 
Although  the  treatment  of  creep  and  consolidation  as  two  distinct  phenomena 
would  seem  appropriate  for  predicting  delayed  settlements  in  frozen  soils, 
experimental  limitations  necessitate  that  these  two  phenomena  be  lumped 
together.  Thus,  frozen  soil  is  treated  as  a  quasi-single-phase  medium  with 
an  empirically  derived  creep  relationship.  Recognizing  that  under  certain 
conditions  consolidation  volume  change  may  have  a  significant  impact  on  the 
delayed  settlement  of  a  footing,  Ladanyi  (1981)  studied  the  frozen  soil’s 
response  to  a  stress  increment  with  the  aid  of  the  Rendulic  plot.  These 
plots  conveniently  enable  any  given  stress  increment  in  a  triaxial  test  to 
be  separated  into  its  hydrostatic  and  deviatoric  components. 

Interaction  between  strength  and  deformation  is  illustrated  by 
considering  a  triaxial  test  on  a  normally  consolidated  sample  of  frozen  soil 
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containing  i hree  phases:  mineral  particles,  ice  and  unfrozen  water.  The 
sample  is  initially  in  equilibrium  at  point  O',  Figure  5.11.  A  significant 
difference  between  frozen  and  unfrozen  soil  is  that  a  total  stress 
increment.  Act ,  (O’A),  can  be  applied  to  a  frozen  soil  exceeding  its  long 
term  strength.  As  the  soil  is  allowed  to  deform  under  closed  system 
conditions,  the  stresses  are  internally  shared  by  the  mineral  particles  and 
the  pore  filling.  The  unfrozen  water  can  only  support  hydrostatic  pressure. 
The  shear  stresses  are  supported  by  the  mineral  phase  and  temporarily  by  the 
ice.  Since  the  ice  phase  supports  a  portion  of  the  shear  stress,  the 
effective  strength  is  mob ilized  only  to  point  B  ’  .  The  hydrostatic  pressure 

generated  by  this  straining,  assumed  to  be  equal  in  the  ice  and  unfrozen 
water,  is  given  by: 


Act 


oct 


-  Act 


oct 


(5.7) 


i  his  is  somewhat  less  than  for  the  same  soil  in  an  unfrozen  state, 
oimilarly,  the  shear  stress,  at  ,  .  acting  in  the  ice  phase  is  given  by: 


At 


oct,  i 


4toc+'4t 


oct 


(5.8) 


where:  AToct  anc^  Axoct  are  *  applied  shear  stress  increment 

and  effective  shear  stress  assumed  by  the  soil,  respectively. 

Maintaining  closed  system  conditions,  the  frozen  soil  will  creep  and  the  ice 
will  gradually  transfer  its  applied  shear  stress  to  the  soil  skeleton.  The 
point  B’  will  move  along  the  effective  stress  path  to  mobilize  full  soil 
strength  at  B.  Creep  brings  the  ice  closer  to  failure  producing  a  loss  of 
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strength  with  fime.  The  failure  surface  shrinks  from  its  initial  position 

:-0  until  failure  occurs  at  point  A  at  t  =  t^.  Since  point  A  lies 

beyond  the  long  term  strength,  creep  failure  will  inevitably  occur.  The 

stresses  at  failure  will  be  Ax'  and  Act1  in  the  soil  skeleton,  au.  in 

oct  oct  i 

the  ice  and  unfrozen  water  and  Ax  .  jn  jce. 

T  *  I 

Opening  the  system  at  point  B’  will  initiate  consolidation  following 
the  stress  path  B'A'.  The  total  deformation  of  the  soil  will  consist  of 
consol idation  and  creep,  occurr i ng  s i mu  I taneous ly.  Full  soil  strength 

mobilization  at  point  A’  will  depend  on  the  rate  of  consolidation  relative 
to  the  steady  state  creep  rate. 

The  foregoing  example  of  a  tr i ax i a  I  specimen  clearly  showed  that  at 
least  for  short  term  conditions,  the  strength  of  the  ice  matrix  can  be 
relied  on  to  support  stress  increments  that  exceed  the  soil’s  long  term 
strength.  However,  the  shape  and  position  of  the  long  term  strength 

envelope  must  be  known  to  clearly  define  the  long  term  response  of  a  frozen 
soil  to  an  arbitrary  stress  increment. 

In  using  the  quasi-single  phase  approach,  it  is  important  that  the 
stress  paths  followed  in  the  I aboratory  to  develop  the  constitutive  equation 
coincide  with  those  expected  under  field  conditions.  The  stress  paths  shown 
in  Figure  5.11  are  for  an  axially  symmetric  stress  state  which  would  be 
suitable  for  calculating  the  time  dependent  response  of  the  frozen  soil 
underneath  a  circular  footing.  Under  conditions  of  plane  strain,  such  as 
the  stress  state  in  the  soil  mass  surrounding  an  underground  opening,  the 
use  of  a  simple  stress  space  in  which  the  coordinates  are  the  vertical  and 
horizontal  principal  stresses  would  be  more  appropriate.  This  arises  from 
the  fact  that  the  intermediate  principal  stress  is  generally  unknown. 
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The  stress  paths  shown  in  the  Rendulic  plot,  Figure  5.11,  have  been 
replotted  in  a  vertical  and  horizontal  stress  space  as  shown  in  Figure  5.12. 
In  this  plot,  the  stress  paths  for  undrained  extension  tests  are  also  shown. 
The  deviatoric  stress  is  given  by  either  the  vertical  or  horizontal  distance 
between  the  point  character i z i ng  the  stress  state  and  the  hydrostatic 
pressure  line,  a  =  o'. 

Although  a  clear  understanding  of  effective  stress  changes  in  a  frozen 
soil  during  shear  has  not  been  developed  to  date,  Ladanyi’s  (1981) 
conceptualization  provides  a  useful  framework  to  study  the  response  of  a 
frozen  soil  to  a  stress  increment.  The  interaction  between  strength  and 
deformation  wi I  I  shed  some  I  ight  on  the  in-situ  deformation  behavior  of  the 
CRREL  room  in  the  Fox  Tunnel. 

Excavation  of  an  underground  opening  will  initiate  deformation  in  the 
soil  mass  surrounding  the  opening.  Considering  the  CRREL  room  at  the  Fox 
Tunnel,  the  ice-rich  silt  comprising  the  roof  will  deform  much  more  readily 
than  the  frozen  dense  gravel  forming  the  walls.  Downward  movement  of  the 
unsupported  ceiling  must  now  be  resisted  by  shear  stresses  that  develop  at 
the  interface  between  the  yielding  and  stationary  soil  mass. 

Before  excavation  of  the  USA  CRREL  room,  the  vertical  stress  along  the  ice- 
rich  silt/gravel  contact  was  equal  to  the  overburden  pressure.  After 
excavation,  the  vertical  stress  along  the  ceiling  of  the  room  is  zero.  This 
stress  release  is  accompanied  by  a  correspond i ng  increase  in  stress  in  the 
adjoining  soil  mass  above  the  frozen  gravel.  Considering  the  two  elements 
shown  in  Figure  5.13(a),  the  stress  in  element  A  will  be  a  small  fraction  of 
what  it  was  before  excavation  commenced.  The  vertical  pressure  in  element  B 
will  have  to  increase  by  the  same  amount  to  maintain  equilibrium.  The 
horizontal  stress  these  two  elements  exert  upon  each  other  must  remain  equal 


. 


. 


* 


105 


and  opposite.  The  induced  shear  stress  at  the  interface  between  elements  A 
and  B  will  cause  a  rotation  of  the  principal  stresses  above  the  tunnel 
supports . 

Although  the  actual  stress  paths  followed  by  the  soil  elements 
surrounding  a  rectangular  opening  are  quite  complicated,  some  insight  into 
the  interaction  between  strength  and  deformation  of  frozen  soil  may  be 
gained  by  studying  the  simple  stress  paths  shown  in  Figure  5.13(b).  The 
stress  release  in  element  A  is  given  by  the  stress  path  O'A  while  an  equal 
but  opposite  stress  increment  is  applied  to  element  B,  ie.  path  O'B.  As 
shown  in  the  figure,  the  stress  release  in  element  A  exceeds  the  long  term 

aH  >  aV  ^a'lure  line.  However,  an  equal  stress  increase  in  element  B  lies 
below  the  long  term  ay  >  aH  failure  line. 

This  simplified  stress  path  description  clearly  shows  the  importance  of 
considering  the  stress  path  in  selecting  an  appropriate  constitutive 
relationship  for  predicting  the  initial  time  dependent  response  of  the 
frozen  soil  surrounding  an  underground  opening.  Since  the  stress  increase 
in  the  soil  above  the  gravel  does  not  exceed  the  soil's  long  term  strength, 
the  creep  behavior  may  well  continue  to  attenuate  with  time.  On  the  other 
hand,  the  soil  above  the  roof  of  the  Tunnel  that  experiences  a  stress 
release  that  exceeds  the  soil's  long  term  strength  will  creep  at  either  a 
steady  state  or  decreasing  rate.  The  creep  rate  will  depend  on  the  relative 
position  of  the  stress  point  A  with  respect  to  the  service  life  strength, 
envelope  for  the  frozen  soil  under  consideration.  Failure  will  inevitably 
occur  for  the  purely  frictional  long  term  strength  shown  in  Figure  5.13(b). 

Evidence  of  the  frozen  silt  approaching  failure  can  be  inferred  from 
the  vertical  closure  rate  versus  time  plot  at  Sta  1+  83.  Installation  of 
the  USBM  bulkhead  at  Sta  1  +  20  increased  the  permafrost  temperature  at  Sta. 
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1+83.  The  warmer  soil  temperatures  resulted  in  a  50%  increase  in  vertical 
closure  rate  at  Sta.  1+83  over  a  period  of  three  months.  A  portion  of  the 
increase  is  due  +o  the  temperature  dependence  of  creep.  However,  the  warmer 
soil  temperatures  will  increase  the  unfrozen  water  content  of  the  frozen 
silt.  Increasing  a  frozen  soils  unfrozen  mo i sture  content  will  shrink  the 
delayed  failure  surface  from  its  original  position  (ie.  instantaneous  or 
short  term  strength)  towards  the  long  term  failure  line.  The  soil 
temperature  becomes  cooler  following  removal  of  the  USBM  bulkhead  at  Sta. 
1+20.  As  the  soil  temperature  decreases,  its  failure  surface  expands 
towards  the  instantaneous  strength.  Thus,  a  numerical  prediction  of  the  i  n- 
s i tu  deformation  behavior  of  the  CRREL  room  would  have  to  account  for 
temperature  dependent  creep  properties  and  a  temperature  and  time  dependent 
failure  criterion  governing  the  strength  of  the  frozen  soil. 

The  vertical  deformation  instrumentation  installed  in  the  CRREL  and 
USBM  rooms  at  the  Fox  Tunnel  was  not  able  to  measure  vertical  deformations 
in  the  gravel  walls  independently  of  the  overall  vertical  closure  of  the 
room.  The  arching  support  provided  by  the  overlying  silt  (due  to  the  stress 
release  above  the  roof)  will  increase  the  vertical  stress  in  the  frozen 
gravel  walls  in  the  vicinity  of  the  room.  The  increase  in  vertical  stress 
coupled  with  a  horizontal  stress  release  brought  about  by  excavating  the 
room  will  significantly  increase  the  shear  stress  in  the  frozen  gravel.  The 
build  up  of  shear  stress  in  the  frozen  gravel  walls  will  lead  to  time 
dependent  deformations  occurring  in  the  gravel  that  may  eventually  lead  to 
de I ayed  fa i I ure. 


Thompson  and  Sayles  (1972)  incorporated  a  yield  criterion  in  their 
finite  element  program  to  account  for  plastic  yield  of  the  material.  In 
their  analysis,  the  effective  stress  as  defined  by  Odqvist  (1966)  was  used 
as  the  yield  criterion.  After  each  creep  increment,  the  effective  stress 
for  each  element  was  checked  to  determine  whether  it  had  exceeded  a  given 
stress  level.  Plastic  deformations  were  allowed  to  take  place  in  all 
elements  that  exceeded  the  given  yield  stress.  It  is  interesting  to  note 
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that  when  the  in  situ  creep  behavior  of  the  frozen  silt  and  gravel  was 
governed  by  primary  creep,  plastic  yielding  was  initiated  in  the  gravel 
walls.  The  elastic-plastic  analysis  carried  out  by  Thompson  and  Sayles 
(1972)  indicate  that  a  catastrophic  failure  would  occur  which  lead  the 
authors  to  conclude  that  the  long  term  strength  of  the  frozen  silt  and 
gravel  was  zero  and  that  the  in  situ  creep  behavior  was  more  accurately 
represented  by  a  steady-state  creep  relationship. 

ine  use  of  effective  stress  f or  a  failure  criterion  is  over-simplified. 
Numerous  studies  have  shown  that  the  strength  of  frozen  soil  is  governed  by 
a  Mohr-Cou I omb  failure  law  where  the  shape  of  the  failure  envelope  depends 
on  soil  type,  density,  ice  saturation,  temperature  and  strain  rate.  The 
strength  of  the  gravel  will  most  certainly  be  governed  by  a  Mohr-Coulomb 
fa i I ure  enve I  ope . 

It  is  possible  that  the  frozen  gravel  was  approaching  the  failure 
envelope  and  experiencing  large  vertical  strains  during  the  period  of  time 
the  soil  temperatures  increased  while  USBM  bulkhead  was  in  place  at  Sta 
1t20.  Thus,  vertical  displacements  in  the  gravel  walls  may  be  superimposed 
on  to  vertical  creep  displacement  in  the  frozen  silt  to  give  the  total  room 
closure  as  recorded  by  Thompson  and  Sayles  (1972).  The  strain  rate  measured 
in  the  laboratory  is  independent  of  the  vertical  deformations,  in  the  gravel 
walls.  This  is  very  likely  why  the  creep  equation  which  gave  the  best  fit 
to  the  observed  room  closured  yielded  a  strain  rate  3.3  times  greater  than 
the  strain  rate  measured  in  the  laboratory. 

5*5  Underground  Circular  Cavities  in  Permafrost 


Major  design  considerations  of  tunnels  in  frozen  ground  are  stand-up 
time,  change  in  diameter  of  the  unlined  opening  with  time  and  the  change  in 
pressure  on  the  tunnel  liner  with  time.  The  first  two  factors  will 
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influence  the  decision  on  whether  or  not  to  support  the  underground  opening 
with  a  tunnel  liner.  The  third  factor  will  provide  guidelines  for 
dimensioning  the  tunnel  liner.  A  small  parametric  study  of  shallow  circular 
tunnels  in  frozen  ground  was  carried  out  to  study  closure  phenomena  and  the 
stress  changes  in  the  frozen  soil  surrounding  the  opening. 

The  first  part  of  this  study  deals  with  the  closure  behavior  of  an 
uni ined  circular  opening  excavated  at  various  depths  in  warm  ice-rich 
permafrost.  The  warm  ice-rich  permafrost  was  chosen  as  the  host  soil  medium 
as  it  would  present  an  upper  bound  for  the  closure  rates.  A  circular 
opening  was  chosen  since  it  is  more  amenable  to  simpler  analysis  under 
certain  conditions. 

The  height  to  diameter  (H/D)  ratios  of  the  circular  openings  varied 
from  1.5  to  4.0.  The  finite  element  mesh  for  H/D=2.0  is  shown  in  Figure 
5.14.  A  homogeneous  soil  profile  was  assumed  for  all  analysis.  Material 
properties  were  assigned  on  the  basis  of  data  reported  in  the  I iterature. 

The  flow  law  for  frozen  soil  used  in  this  study  was: 

I  =  2.0  x  ,0-V  (5.9) 

This  flow  law  is  identical  to  the  flow  law  for  ice  at  -2°C  (Morgenstern 
et  al.,  1980)  and  represents  a  valid  upper  bound  to  the  in-situ  creep 
behavior  of  fine-grained,  ice-rich  frozen  soils  of  the  Mackenzie  River 
Val ley,  N.W.T. 

A  plot  of  the  closure  velocity  versus  depth  of  overburden  for  the 
crown,  springline  and  invert  of  the  tunnel  is  shown  in  Figure  5.15.  The 
depth  of  overburden  in  this  figure  is  expressed  non-d imens iona I  I y  as  the  H/D 
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ratio,  ^or  al I  three  points  considered,  the  closure  rate  increases  with  the 
depth  of  overburden.  Also,  the  downward  movement  of  the  crown  is  greater 
than  the  upward  movement  of  the  invert  creating  a  sag  in  the  roof.  Thus, 
for  shallow  tunnels,  the  initial  circular  cross  section  undergoes  a 
transition  to  an  oval  shape  as  the  frozen  soil  creeps  inward.  At  greater 
depths,  this  transition  will  be  less  pronounced  because  the  relative  change 
in  overburden  stress  from  crown  to  invert  wi  I  I  be  smal  ler.  Based  on  these 
preliminary  results,  it  would  appear  that  a  circular  opening  located  9  to  10 
diameters  below  the  surface  would  behave  axisymmetr ica I  I y .  This  is  somewhat 
deeper  than  a  circular  opening  in  an  elastic  soil  or  rock. 

Velocity  vectors  depicting  magnitude  and  direction  throughout  the 
region  analyzed  for  H/D=2.0  and  H/D=3.0  are  shownin  Figures  5.16  and  5.17, 
respectively.  The  relative  magnitude  in  these  figures  is  given  by  the 
length  of  the  line  segment.  For  both  tunnel  depths,  the  frozen  soil  flows 
in  towards  the  opening.  The  transformed  cross  section  is  shown  on  a  grossly 
exaggerated  scale  by  the  dashed  line  joining  the  velocity  vectors  around  the 
periphery  of  the  opening. 

Contours  of  effective  strain  rate  for  H/D=2.0  and  H/D=3.0  are  shown  in 
Figures  5.13  and  5.19,  respectively.  Although  the  magnitude  of  individual 
strain  rate  components  are  not  given  by  the  effective  strain  rate,  the 
contours  serve  the  purpose  of  indicating  where  the  high  strain  rate 
gradients  exist.  The  greatest  effective  strain  rate  gradients  directly 
above  the  crown  of  the  opening  and  diminish  as  one  sweeps  around  from  the 
crown  to  the  invert.  Higher  gradients  exist  for  the  tunnel  located  at 
H/D=3.0  (see  Figure  5.19)  indicating  •"'hat  the  disturbing  influence  of  the 
ground  surfaces  decreases  with  deeper  tunnels.  Both  Figures  5.18  and  5.19 
show  the  high  strain  rate  gradients  extend  only  1.5  to  2  tunnel  diameters 
above  the  crown  of  the  opening.  This  zone  of  influence  will  decrease  with 
deeper  tunnels. 
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11"  was  c  I  ear  from  1"he  foregoing  analyses  that  shal  low  circular  cavities 
located  in  warm  ice-rich  permafrost  would  become  unstable  unless  some 
measures  were  undertaken  to  arrest  the  deformations.  The  second  part  of 
this  study  dealt  with  the  closure  of  a  I ined  circular  opening  in  permafrost. 
Only  one  opening  located  2  diameters  below  the  surface  was  analyzed.  The 
same  finite  element  mesh  as  shown  in  Figure  5.14  was  used  except  that  a 
narrow  ring  of  elements  was  added  to  the  periphery  of  the  opening  to 

simulate  a  100  mm  thick  concrete  liner.  The  geometry  is  shown  in  Figure 
5.20. 


This  analysis  was  carried  out  by  calculating  the  elastic  gravity  stress 
field  surrounding  the  circular  opening  without  the  I  iner  elements  present. 
After  this  step,  the  elements  representing  the  tunnel  liner  were  added  to 
the  geometry  and  the  creep  solution  was  initiated.  It  was  felt  that  this 
sequence  would  be  more  representat i ve  of  actual  construction  operations. 

The  velocity  vector  diagram  for  the  I ined  tunnel  placed  2  diameters 
below  the  ground  surface  is  shown  in  Figure  5.21.  As  shown  in  the  figure, 
the  tunnel  liner  plays  a  significant  role  in  altering  the  velocity 
distribution  in  the  frozen  soil  surrounding  the  opening.  The  overall 
movement  of  the  soil  mass  and  tunnel  opening  is  downward  due  to  gravity 
stress.  In  the  vicinity  of  the  tunnel,  the  creeping  frozen  soil  is 
deflected  around  the  tunnel  since  its  stiffness  is  much  greater  than  that  of 
the  frozen  soil.  The  velocity  profile  is  a  function  of  the  stiffness 
properties  of  the  tunnel  liner.  The  tunnel  liner  also  serves  to  drastically 
reduce  the  downward  velocity  of  the  crown  of  the  tunnel.  In  this  case,  the 
crown  velocity  is  reduced  by  3  orders  of  magnitude  for  a  tunnel  at  the  same 
depth  with  no  liner  present. 
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FIGURE  5.1 


SITE  PLAN  OF  FOX  TUNNEL  SITE, 
ALASKA  (Sellman,  1967) 
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FIGURE  5.2  PLAN  FOF  FOX  TUNNEL  SFIOWING  SECOND  STAGE  OF 

EXCAVATION  (Pettibone,  1973) 
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FIGURE  5.3 


GEOLOGIC  CROSS  SECTION,  FOX  TUNNEL  SITE 
(Pettibone,  1973) 
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FIGURE  5.4 


CLOSURE  MEASUREMENTS  OF  MAIN  TUNNEL 
FOLLOWING  1965  CONSTRUCTION  SEASON 
(Swinzow,  1970) 
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FIGURE  5.5 


CLOSURE  MEASUREMENTS  OF  MAIN  TUNNEL 
FOLLOWING  1965  CONSTRUCTION  SEASON 
(McAnerney,  1968) 
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FIGURE  5.6 


VERTICAL  CLOSURE  OF  USBM  ROOM  IN  THE  FOX  TUNNEL 
(Pettibone,  1973) 
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FIGURE  5.7  VERTICAL  CLOSURE  OF  USA  CRREL  ROOM  IN  THE  FOX  TUNNEL 
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FIGURE  5.8  HORIZONTAL  CLOSURE  OF  USA  CRREL  ROOM  IN  THE 
FOX  TUNNEL,  ALASKA 
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(a)  VERTICAL  CLOSURE  RATE 
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FIGURE  5.9  VERTICAL  AND  HORIZONTAL  CLOSURE  RATE  OF  USA  CRREL  ROOM 

IN  THE  FOX  TUNNEL,  ALASKA 
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FIGURE  5.10  FINITE  ELEMENT  DISCRETIZATION  OF  USA  CRREL  ROOM 

IN  THE  FOX  TUNNEL 
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FIGURE  5.11 


RENDULIC  PLOT  OF  NORMALLY  CONSOLIDATED 
PLASTIC  FROZEN  SOIL  STRESSED  BEYOND  ITS 
LONG  TERM  STRENGTH  (Ladanyi,  1981) 


FAILURE  LINES 

aV  >aH 


FIGURE  5.12 


STRESS  PATHS  FOR  NORMALLY  CONSOLIDATED 
PLASTIC  FROZEN  SOIL  STRESSED  BEYOND  ITS 
LONG  TERM  STRENGTH 
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FIGURE  5.13  STRESS  PATHS  FOR  SOIL  ELEMENTS  ABOVE  AN 

UNDERGROUND  OPENING 
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FIGURE  5.14 


FINITE  ELEMENT  MESH  FOR  TUNNEL  AT  H/D  =  2.0 
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FIGURE  5.15 


TUNNEL  CLOSURE  VELOCITY  VERSUS  DEPTH  OF  OVERBURDEN 
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FIGURE  5.16 


VELOCITY  VECTOR  PLOT  FOR  TUNNEL  AT  H/D  =  2.0 
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FIGURE  5.17 


VELOCITY  VECTOR  PLOT  FOR  TUNNEL  AT  H/D  =  3.0 
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FIGURE  5.18  EFFECTIVE  STRAIN  RATE  CONTOURS  FOR  TUNNEL  AT  H/D  =  2.0 
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FIGURE  5.19 


EFFECTIVE  STRAIN  RATE  CONTOURS  FOR  TUNNEL  AT  H/D  =  3.0 
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FIGURE  5.20 


FINITE  ELEMENT  MESH  FOR  LINED  TUNNEL  AT  H/D  =  2.0 
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10~3  mm/yr 


FIGURE  5.21 


VELOCITY  VECTOR  PLOT  FOR  LINED  TUNNEL  AT  H/D  =  2.0 


CHAPTER  VI 


CONCLUD I NG  REMARKS 


6.1  GENERAL 


An  analytical  study  of  the  in-situ  creep  behavior  in  ice-rich  frozen 
soil  has  been  reported  herein.  A  review  of  the  deformation  behavior  of 
frozen  soil  has  lead  to  the  conclusion  that  the  steady  state  creep  of  ice- 
rich  frozen  soil  is  adequately  described  by  an  empirically  derived  simple 
power  law  which  relates  strain  rate  to  stress.  The  flow  law  for  steady 
state  creep  in  polycrystalline  ice  (Morgenstern  et  al.,  1980)  represents  a 
valid  upper  bound  for  the  constitutive  behavior  of  frozen  soil  under 
equivalent  loading  conditions.  The  importance  of  testing  natural  permafrost 
samples  with  representat i ve  segregated  ice  features  was  also  noted. 

The  incremental  initial  strain  finite  element  procedure  was  applied  to 
two  case  histories  of  naturally  occurring  creep  in  ice-rich  permafrost.  The 
in-situ  creep  behavior  of  a  relatively  steep  slope  on  the  left  bank  of  the 
Great  Bear  River  at  the  proposed  Arctic  Gas  crossing  was  studied  by  Savigny 
(1980).  Naturally  occurring  creep  in  ice-rich  Fairbanks  silt  at  the  Fox 
Tunnel  near  Fairbanks,  Alaska  was  studied  by  Swinzow  (1970),  Thompson  and 
Sayles  (1972)  and  Pettibone  (1973).  The  following  sections  summarize  the 
major  conclusions  of  the  analytical  creep  studies  at  these  two  sites. 
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6.2  Cone  I  us  ions 


6.2.1  Great  Bear  River 

The  finite  element  creep  analysis  of  the  left  bank  of  the  Great  Bear 
River  at  the  proposed  Arctic  Gas  crossing  has  shown  that  the  in-situ  creep 
behavior  of  the  g I ac io I acustr i ne  clay  can  be  accurately  represented  by  the 
simple  power  law.  The  flow  law  which  gave  the  best  fit  between  observed  and 
predicted  velocities  had  an  exponent  of  3.0  and  a  coefficient  of 

3.33  x  10  9.  This  strain  rate  is  six  times  slower  than  the  value  for 
polycrystalline  ice  at  -2°C  (Morgenstern  et  al.,  1980). 

The  exact  form  of  the  constitutive  relationship  for  the  g I ac iode I ta ic 
sand  remains  inconclusive.  In  the  first  creep  analysis,  the  creep  behavior 
of  the  sand  was  assumed  to  be  given  by  the  simple  power  law  with  the  same 
parameters  as  those  in  the  g I ac io I acustr i ne  clay.  The  velocities  predicted 
by  the  simple  Dower  law  were  much  greater  than  the  observed  values.  A 
second  analysis  was  carried  out  in  which  the  sand  was  treated  as  an  elastic 
material.  This  proved  to  be  much  too  serious  a  restriction  on  the  creep 
movements  in  the  sand  and  the  g I ac io I acustr i ne  clay  where  it  was  overlain  by 
sand.  This  analysis  did  predict  the  build  up  of  horizontal  tensile  stress 

in  the  sand  much  greater  than  the  sand  would  be  able  to  sustain  under  in- 
situ  conditions. 

The  creep  movements  recorded  at  the  GB1 A  and  GB3  inclinometer  locations 
at  the  Great  Bear  River  showed  that  more  movement  was  associated  with  large 
ice  lenses  when  they  are  widely  separated  (Savigny,  1980).  Numerical 
modelling  of  permafrost  as  composite  soil- ice  medium  would  require  an 
accurate  mapping  of  the  lateral  continuity  of  these  large  ice  lenses. 
Statistical  analysis  of  the  data  derived  from  the  MVPRL  ice  variability 
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test  sites  at  Norman  Wells  N.W.T.  and  Landing  Lake  N.W.T.  suggest  that  for  a 
uniform  s trat i grap h i c  unit,  the  overal I  ground  ice  content  of  permafrost  can 
be  described  as  a  normally  distributed  variable.  Thus,  for  practical 
purposes,  the  assumption  of  homogeneity  is  va! id  for  predicting  overal I 
t i me— dep en den t  ground  movements  in  ice— rich  frozen  soi I. 

A  general  conclusion  derived  from  the  study  of  the  left  bank  of  the 
Great  Bank  River  at  the  proposed  Arctic  Gas  crossing  that  the  time-dependent 
deformation  behavior  of  homogeneous  slopes  consisting  of  ice-rich  fine¬ 
grained  permafrost  can  be  model  led  by  a  simple  power  law.  Using  the  creep 
parameters  for  polycrystal  I ine  ice  wi I  I  provide  an  upper  bound  to  the  creep 
velocities. 

6.2.2  Fox  Tunnel 

The  in-situ  deformation  studies  of  the  Fox  Tunnel  reported  by  Swinzow 
(1970),  Thompson  and  Sayles  (1972),  and  Pettibone  (1973)  showed  that 
underground  openings  excavated  in  warm  ice-rich  permafrost  can  become 
unstable  when  the  ambient  air  temperature  is  near  the  soil  temperature.  In 
each  case,  primary  creep  behavior  contributed  substantially  to  the  overall 
room  closure. 

The  closure  data  of  the  USA  CRREL  room  recorded  over  a  3  year  period 
clearly  showed  attenuating  creep  for  well  over  one  year.  The  closure  versus 
time  data  for  the  USA  CRREL  room  and  the  USBM  room  in  the  Fox  Tunnel  both 
have  the  same  general  form  for  the  first  3  years  of  operation.  Thus, 
Thompson  and  Sayles  (1972)  were  premature  in  concluding  that  steady  state 
creep  dominated  the  closure  of  the  USA  CRREL  room  during  the  first  year  of 
operation. 
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A  finite  element  creep  analysis  of  the  USA  CRREL  room  showed  that  the 
simple  power  law  for  polycrystalline  ice  does  not  represent  a  valid  upper 
bound  for  the  in-situ  creep  behavior  of  the  frozen  Fairbanks  silt.  However, 
analyzing  the  stress  changes  above  a  rectangular  underground  opening  in 
terms  of  arching  theory  clearly  shows  that  the  stress  paths  followed  by  soil 
elements  above  the  roof  of  the  opening  are  entirely  different  from  that 
normally  used  in  constant  stress  or  constant  strain  rate  creep  tests.  The 
strain  rate  predicted  by  finite  element  computations  was  3.3  times  faster 
than  the  laboratory  measured  value  (Thompson  and  Sayles,  1972)  because  the 
vertical  room  closure  measurements  could  not  distinguish  between  vertical 
compression  in  the  gravel  walls  and  downward  creep  movement  of  the  overlying 
ice-rich  silt.  For  this  particular  case  it  was  not  possible  to  isolate  in- 
situ  creep  deformations  by  measuring  the  closure  between  opposing  points  in 
the  roof,  floor  and  walls. 

6*3  Recommendations  for  Further  Study 

The  analytical  study  of  naturally  occurring  creep  in  a  slope  and 
underground  cavity  located  in  ice-rich  permafrost  has  shown  that  the  simple 
power  law  can  be  used  to  predict  in-situ  deformations  provided  that  the 
stress  does  not  approach  failure. 

The  largest  discrepancies  between  predicted  and  observed  behavior 
appear  to  coincide  with  regions  where  the  stresses  experienced  by  the  soil 
mass  are  unlike  the  stress  paths  used  in  laboratory  constant  stress  or 
constant  strain  rate  creep  tests.  This  was  illustrated  at  the  GB2 
inclinometer  location  at  the  Great  Bear  River.  The  stresses  at  this 
location  are  altered  from  the  simple  shear  condition  because  the  downslope 
creep  movement  in  the  g I aciolacustr ine  clay  is  impeded  by  the  overlying  sand 
just  upslope  of  this  location.  Downslope  of  the  GB2  inclinometer,  the 
glaciolacustr ine  clay  creeps  at  a  faster  rate  since  there  is  no  overlying 
sand  to  restrict  movement.  Thus,  the  clay  in  the  vicinity  of  the  GB2 
inclinometer  would  experience  an  unloading  in  the  horizontal  direction. 
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The  frozen  silt  above  the  roof  in  the  USA  CRREL  room  in  the  Fox  Tunnel 
experiences  an  unloading  due  to  excavation  of  the  room.  Arching  behavior 
will  transfer  the  vertical  stress  release  in  the  roof  to  the  frozen  silt 
above  the  walls  and  to  the  frozen  gravel  walls.  Thus,  while  the  soil  which 
experiences  greater  compressive  stress  may  still  lie  below  the  compression 
failure  envelope,  the  soil  above  the  roof  of  the  tunnel  may  experience 
tensile  failure  at  some  time  after  the  room  is  excavated.  The  stress  state 
can  become  guife  complicated  whenever  two  or  more  frozen  materials  are 
present  with  different  creep  properties.  Thus,  it  is  important  to  obtain  an 

accurate  assessment  of  the  creep  behavior  of  the  materials  relative  to  each 
other. 

Future  numerical  modelling  of  ice-rich  permafrost  under  complex  stress 
conditions  should  be  directed  towards  studying  the  stress  changes 
experienced  by  the  soil  as  deformations  develop  and  to  incorporate  the 
concept  of  a  limiting  long  term  strength  into  the  analysis.  Toward  this 
end,  laboratory  studies  should  be  carried  out  to  study  the  creep  behavior  of 
frozen  soils  under  stress  conditions  that  simulate  unloading.  Also,  future 
field  programs  should  have  carefully  designed  and  installed  instrumentation 
to  ensure  that  the  in-situ  creep  deformations  can  be  isolated  from  other 
deformation  processes  which  may  be  occurring  simultaneously. 
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APPENDIX  A 

FORMULATION  OF  FINITE  ELEMENT  EQUATIONS  FOR 
INCREMENTAL  INITIAL  STRAIN  PROCEDURE 


This  appendix  presents  a  brief  review  of  how  the  incremental  theory  of 
creep  is  introduced  into  the  finite  element  method  using  the  initial  strain 
procedure. 

In  the  incremental  theory  of  creep,  the  increment  of  creep  strain  is 
treated  as  an  initial  strain  for  any  one  time  interval  and  is  assumed  to  be 
constant  for  that  time  interval. 

The  increment  of  total  strain,  AeT,  for  any  time  interval  can  be 
assumed  to  consist  of  an  increment  of  elastic  strain,  Ae^  ancj  an  jncrement 
increment  of  creep  strain,  Ae^.  This  can  be  expressed  as: 

{  Ae^"}  =  {  Ae  +  (AeS  ( /\]  ) 


Solving  equation  A1  for  the  elastic  strain  gives: 


{Ae^}  =  (Ae^"}  -  {Ae^} 


( A2 ) 


and  applying  Hooke’s  law: 


{Aa}  =  [ D 1 ( AeT  -  AeC} 


(A3) 


where: 


[ D 1  =  matrix  of  elastic  constants 


. 


■ 
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For  rhe  linear  elastic  case,  an  increment  of  strain  energy  density  of  a  body 
is  g i ven  by : 


AU(x,y,z)  =  1  {AeE}T{Aa} 
2 


( A4 ) 


Substitution  of  equation  A2  and  A3  into  A8  gives: 


AU(x,y,z)  =  J_  {AeT  -  A£C}T[DHAeT  -  AeC} 
2 


( A5 ) 


An  increment  of  potential  energy  of  a  body  is  expressed  mathematically  as: 


Air  =  /yAU(x,y,  z)  dV  — /y(X  Au  +  YAv  +  ZAw)dV 


—  Y  ( T  Au  +  I  AV  t  T  Aw)  dS-p 
S  x  y  z  T 


( A6 ) 


where  V  represents  the  volume  of  the  body  and  ST  is  +he  surface  on  which 

surface  tractions  are  prescribed.  The  last  two  integrals  in  equation  A6 

represent  the  work  done  by  the  external  forces;  i  .e.  the  body  forces  X,  Y, 

Z,  and  the  surface  tractions,  T  ,  T  ,  and  T  .  The  bar  denotes  quantities 

x  y  z 

that  are  prescribed.  Substituting  equation  A5  into  A6  yields: 


Att  -  1  /  {AeT}T[D]{AeT}dV  -  1  /  {  AeT}  T  [  D]  { AeS  dV 
p  2  V  2  V 


Q  j  -p 

-  4  -MAe  }  t  D  ]  {  Ae  }  dV  +  1  /  {Ae'VtDHAeS  dV 
2  V  2  V 


( A7 ) 


-  1  /  {  A4)  ( x,  y ) }  T{  F}  dV  -  /  {  AiJj  ( x,  y )  T{  T}  dS 

2  V  ST  T 


where:  {  A^(x,  y ,  z) }  T=  [  Au  Ay  Az  J 


Using  standard  finite  element  calculations,  one  can  write: 

{Ae}  =  [ B I  {Aq}  (A 8) 

{ At ( x, y , z) }  =  INI  {Aq}  (A9) 

where  [B1  =  matrix  relating  element  strains  to  the  values  of  the 

displacement  at  its  nodes. 


[  N I 


coefficient  matrix  for  displacement  interpolation  model. 


'  X 


Substituting  equation  A8  and  A9  into  equation  A7  yields: 


Ait  -If  {  Aq}  T  t  B  I T  [  D  ]  t  B  ]  {  Aq}  dV  -  /  { A q}  T  [  8  1 T  [  D  ]  {  AeS  dV 

p  2  V  y 


+1  f  {AeC}  1  [ 0 ] { Ae C} d V  -  /  { Aq} T [ N ) T{ F} dV 
2  V  y 


-  /$  { Aq} T [ N ] T{ T} dS^ 


( A 1  0  ) 


The  Theorem  of  Minimum  Potential  Energy  states: 

"Of  all  possible  displacement  configurations  a  body  can  assume  which 
s  at i s  f y  comp  at ab i I i i y  on  the  con  s tr a i nts  or  k i nemat i c  bou  n  dar y 
conditions,  the  configuration  satisfying  equilibrium  makes  the 
potential  energy  assume  a  minimum  value". 

Mathematically,  this  is  expressed  as: 


6  (  Att  )  -  0 
P 


(A1  1  ) 


Taking  the  first  variation  of  equation  A10  with  respect  to  the  incremental 
displacements,  {Aq},  and  setting  it  equal  to  zero  gives: 


{6  (Aq)}T[/y[B]T[D]  [  B  ]  d  V  { A  q }  -  /  y  [  B  1  T  [  D 1  (  Ae  C}  dV 

-  /  [  N 1  T(  F>  dV  -  /  [  N1  T(  T>  dS  = 

V  S 


(A12) 


0 


.  2 


. 


Since  the  variation  of  the  displacements  are  arbitrary,  the  quantity 
enclosed  in  the  brackets  must  vanish.  The  variation  of  the  third  integral 
in  equation  A1 0  is  zero  because  it  is  not  a  function  of  the  nodal 
displacements.  The  set  of  equilibrium  equations  for  each  time  interval  is: 


[K1  { Aq}  =  {P}  +  (AF  } 

s  s  s  c  s 


( A 1  3  ) 


where: 


M  el  M  j 

IKJ  =  2  [Ki  =2.  (/  IB1 1 [ D 1 [ B ] d V ) 
s  m=  1  m  rn=  1  V  m 


(P> 

s 


M  e  I  IM 

=  V  P>  =  2  ,(/..(  N )  (  F>  dV  +  f0  [  N 1  (  T>  dS_ ) 
m=1  m  m=1  V  T  m 


M  e  I  VI 

(AF  }  =  2  { A F  }  =  2 

c  s  m=  1  c  m  m=  1 


(/v[B]T[DHAeC}dV) 


m 


The  summation  sign  is  carried  over  the  total  number  of  elements,  M,  to 
obtain  the  equilibrium  equations  for  the  entire  assemblage. 


The  vector  {Fc>s  represents  the  creep  strain  nodal  load  vector. 
Thus,  it  can  be  seen  that  for  each  interval  of  time,  creep  strains  are 
allowed  to  take  place.  Equivalent  nodal  forces  are  then  calculated  which 
would  be  necessary  to  cause  elastic  strains  of  the  same  magnitude.  These 
fictitious  nodal  forces  are  then  added  to  the  load  vector  and  new  set  of 
nodal  displacements  are  calculated  during  that  time  interval  with  equation 
A1  3. 
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CALCULATION  OF  FICTITIOUS  NODAL  FORCES 


This  appendix  presents  the  derivation  of  the  fictitious  nodal  creep 
forces  used  in  the  initial  strain  approach  to  solving  creep  problems. 

In  each  time  interval,  the  increment  of  creep  strain  is  given  by  the 
stress-stra i n-rate  law,  i.e.: 


Ae  .  .  -  3  ( Ae  /a  )  s 
1 J  2  e  e  i  j 


(B1  ) 


Using  Hooke's  law,  the  increments  of  creep  strain  are  converted  to 
incremental  creep  stresses: 


{  qC}  =  [  D 1  { e  C} 


(B2) 


The  constant  stress  state  for  the  constant  strain  triangular  finite  element 
is  represented  by  Figure  8-1 (a). 


The  stresses  can  be  replaced  by  a  set  of  statically  equivalent  forces 
acting  at  the  node  points  of  the  constant  strain  triangle.  Elements  with 
higher  order  displacement  interpolation  functions  require  that  the  nodal 
forces  be  work-equivalent  loads.  However,  for  the  constant  strain  triangle, 
work-equivalent  and  stat i ca I  I y-equ i va I ent  loads  are  equal. 

The  incremental  stresses  are  related  to  the  nodal  forces  for  each 
e I ement  as  fo I  lows  : 


( B3 ) 


■  X 


1 


The  fictitious  nodal  load  vector  is  computed  for  each  element  in  turn  and 
then  summed  to  obtain  the  incremental  load  vector  for  the  entire  assemblage. 
The  incremental  load  vector  is  then  used  to  calculate  a  new  set  of 
incremental  displacements  for  the  current  time  step. 
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a)  CONSTANT  STRESS  DISTRIBUTION 


Fyk 
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b)  STATICALLY-EQUIVALENT  NODAL  FORCES 


FIGURE  B.  1 


STATICALLY  EQUIVALENT  NODAL  FORCES  FOR 
CONSTANT  STRAIN  TRIANGLE 
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INPUT  DATA  FORMAT  FOR  PROGRAM:  CREEP 


if.  -M  5fc  : 


CARD  TYPE  NO.  1  -  TITLE  OF  THE 


PROBLEI 


[ Format ( 20A4 ) ] 


Columns  Variable 
*  *  *  *  *  *  * 


Entry 

*  5)C  5|e  5^ 


1"80  TITLE  Enter  the  heading  information 

that  is  to  be  printed 
with  the  output 


-any  alphanumeric  characters  can  be  entered 


T ^ 


CARD  TYPE  NO.  2  -  PROBLEM  SOLUTION  CONTROL  [ Format ( 2A8 ) ] 


Co  1 umns 

=tc >f:  jf:  :* 

1-8 


9-16 


Var i able 

Entry 

*  *  *  *  * 

CODE  1 

START 

-  For  an  initial  run  of 
the  program,  all  data 
is  input  on  cards 

RSTART 

-  For  a  continuation  of 
a  previous  problem, 
that  was  terminated 
ear  1  y 

C0DE2 

RUN  - 

For  a  normal  run  of  the 

program 

CONTIN  -  Write  the  results  of 
the  final  time  step 
on  to  tape  for 
restarting  the 
problem  at  a  future 
date 


-if  CODE  1  =  RSTART ,  card  sets  1,  2,  3,  11,  12,  14  and  15 
are  required 


-to  restart  a  previous  problem,  the  results  of  the  last  increment 
of  the  terminated  run  are  read  from  unit  ITAPE 


. 
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:  3*c  afc  s+c  5*: 


CARD  TYPE  NO.  3  -  INPUT  AND  OUTPUT  DESIGNATION  [ Forma t ( 2 1 6 ) ] 


:  sfc 


Columns  Variable 
*  *  *  *  *  *  *  *  *  *  *  *  *  *  * 


Entry 

afc  afc  sfc 


1-6  ITAPE 

7-12  JTAPE 


Unit  number  from  which  the 
results  of  a  terminated  run 
will  be  read 

Unit  number  onto  which 
the  results  for  the  final 
time  step  will  be  written 
for  a  run  that  will  be 
terminated  early 


-for  an  initial  run  of  the  program,  all  the  input  data  is 
assumed  to  be  read  from  cards 


-to  continue  a  problem  at  a  future  date,  a  unit  number 
for  JTAPE  must  be  specified;  otherwise  solution  of  the 
problem  will  terminate 


•  N  ■  M 
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*  * 


.CARD  TYPE  NO.  4  MASTER  CONTROL  CARD  [ Format ( 7 1 6 ) ] 


Co  1 umns 

Var i ab 1e 

Entry 

*  *  *  *  *  *  * 

3)C  5^  3tc  Jfc 

1-6 

NN 

Total  number  of  nodes 
(maximum  of  200) 

7-12 

NE 

Total  number  of  elements 
(maximum  of  300) 

13-18 

NUE 

Total  number  of  materials 
for  the  elastic  solution 
(maximum  of  20) 

19-24 

NUC 

Total  number  of  materials 
for  the  creep  solution 
(maxi mum  of  20 ) 

25-30 

NGI 

Gravity  loading  condition 
NGI=0  no  gravity  load 

NGI = 1  gravi ty  load 
cons i dered 

31-36 

NLC 

External  load  condition 

NLC=0  no  external  loads 
NLC=1  external  loads 
cons i dered 

37-42 

KV 

Creep  condition 

KV=0  only  an  elastic 

solution  is  performed 
K V  = 1  creep  solution  is 
required 

■ 


I 
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CARD  TYPE  NO.  5  -  NODAL  POINT  DATA  [ Forma t ( 3 1 5 , 2F 1 0 . 3 ) ] 


I***:*:*::*:*:**:*:*: 


:  ■*  *.  X  *  *  X  * 


Columns  Variable 
*  *  *  *  *  *  *  *  *  *  *  *  *  *  * 


Entry 

*  *  5+C  *  * 


1-5 

6-10 


11-15 


16-25 

26-35 


I 

ND (1,1) 


ND( I ,2) 


X  { I  ) 
Y  ( I  ) 


Node 

X-direction  boundary 
condition  code  for  node  I 
ND ( I , 1 ) =0  f i xed  i n 

X-di rect ion 
ND (  I  ,  1  )  =  1  free 

Y-direction  boundary 
condition  code  for  node  I 
ND ( I , 2 ) =0  f i xed  i n 

Y-direction 
ND ( I , 2 ) = 1  free 

X  coordinate  for  node  number  I 

Y  coordinate  for  node  number  I 


-one  node  per  card 
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*********;*:: 

CARD  TYPE  NO.  6  - 


ELASTIC  MATERIAL  PROPERTIES 
[ Format (F12.0.2F6.3)  ] 


Co  1 umns 
*  *  *  *  *  *  * 

Variable 

st:***:*:;*;:*:* 

Entry 

***** 

1-12 

El  ( I ) 

Young's  Modulus 
number  I 

of 

mater i a  1 

13-18 

U I  (  I ) 

Poi sson' s  Ratio 
number  I 

of 

mater i a  1 

19-24 

UNI  ( I ) 

Uni t  weight  of  i 
number  I 

mater i a  1 

-there  should  be  one  card  for  each  material  up  to  a  maximum  of 


CARD  TYPE  NO.  7  -  ELEMENT  INCIDENCES  [ Forma t ( 6X , 5 16 ) ] 


*********** 


Columns 

Variable 

Entry 

******* 

******** 

***** 

7-12 

N I 

Corner  node  number  I 

13-18 

NJ 

Corner  node  number  J 

19-24 

NK 

Corner  node  number  K 

25-30 

MET 

Elastic  material  number 
indicator  for  each  element 

31-36 

MMET ( I ) 

Creep  material  indicator  for 
each  element 

*****************************************;,.; 


-one  element  per  card 

-node  numbers  I,  J,  K  are  entered  in  a  counter-clockwise 
direction:  the  starting  node  is  arbitrary. 


■ 


• 
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CARD  TYPE  NO.  8  -  EXTERNAL  LOADING  [ Forma t ( 1 6 ) ] 


:  *  *  *  *  * 


Columns  Variable 

*******  *  *  ****** 


Entry 
*  *  *  *  * 


1-6  NNL 


Number  of  noads  loaded 


:  *  *  *  * 


if  NLC-0,  card  type  no.  7  as  well  as  card  type  no.  8  must 
not  appear . 


'  3iC  34"  34;  34;  34c  ^4c  *A*  ^  *X* 

CARD  TYPE  NO.  9  -  EXTERNAL  LOADS  [ Forma t ( I  6 , 2F 1 2 . 3 ) ] 

:  * 


Columns  Variable 

*******  ******** 


Entry 
*  *  *  *  * 


1-6  JNU 

7-18  FF  (  1  ) 

19-30  FF ( 2 ) 


Node  number  of  loaded  node 
Load  in  X-direction 
Load  in  Y-direction 


-one  node  per  card 

-this  card  set  must  not  appear  if  NLC=0 
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*  *  : 


CARD  TYPE  NO.  10  -  NODAL  PRINT  OUT  CONTROL  [ Forma t ( 1 2 1 6 ) ] 


Co  1 umns 


Variable 

**:+:***** 


Entry 

#  sf:  *  5+:  jf: 


1-6 


7-12 


13-18 


67-72 


NNPR ( I ) 


NNPR ( 1+1 ) 


NNPR ( 1+2) 


NNPR ( I +5 ) 


Nodal  print  out  control 
indicator  for  node  number  I 
NNPR(I)=  0  print  out  and 


NNPR (  I  )  =  -  1 
NNPR (  I  )=  K 


no  rotation 
do  not  print 
print  out  and 
rotate  an  angle 
of  K  degrees. 

(K=an  integer) 


Nodal  print  out  control 
indicator  for  node  number  1+1 

Nodal  print  out  control 
indicator  for  node  number  1+2 


Nodal  print  out  control 
indicator  for  node  number  1+5 


-twelve  nodes  per  card 


* 
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CARD  TYPE  NO.  11  -  ELEMENT  PRINT  OUT  CONTROL  [(Format  (1216)] 


^  ■ 


Columns  Variable 
*  *  *  *  *  *  *  *  *  *  *  *  *  *  * 


Entry 


1-6  NEPR(I) 


7-12  NEPRU  +  1) 

13-18  NEPR ( 1+2 ) 

67-72  NEPR ( 1+5 ) 


Element  print 
indicator  for 
number  I 
NEPR ( I  )  =  0 

NEPR ( I ) =  -  1 
NEPR (  I  )=  K 


Element  print 
indicator  for 
number  1+1 

Element  print 
indicator  for 
number  1+2 


Element  print 
indicator  for 
number  1+5 


out  control 
e 1 ement 

Print  out  and 
no  rotation 
do  not  print 
print  out  and 
rotate  an  angle 
of  K  degrees 
( k=an  i nteger ) 

out  control 
e 1 ement 


out  control 
e 1 ement 


out  control 
e 1 ement 


-twelve  elements  per  card 

-if  KV  =  0  (i.e.,  an  elastic  solution  only)  no  more  data  is 
requi red 


-if  a  creep  solution  is  desired,  the  following  card  sets  must 
be  i nc 1 uded . 
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CARD  TYPE  NO.  12  -  CREEP  SOLUION  CONTROL  CARD 

[Format(F12. 0,316) ] 


Columns 

Variable 

******** 

Entry 

1-12 

TDM 

Maximum  time  allowed  for  creep 
solution 

13-18 

MNI 

Maximum  number  of  creep 
i ncrements 

19-24 

NCPR 

Creep  solution  print  control 
NCPR=0  print  all  output 
NCPR=1  print  selective 
output 

25-30 

NSPO 

Total  number  ofcreep 

increments  to  be  printed 


CARD  TYPE  NO.  13  -  CREEP  SOLUTION  PRINT  CONTROL  [Format  (1216)] 


Co  1 umns 

Var i ab le 

Entry 

******* 

******** 

***** 

1-6 

I  PR  I N  (  1 ) 

The  increment  number  for  which 

the  first  set  of  creep  results 
are  to  be  printed 

I  PR  I N ( I  MAX ) 


The  increment  number  for  which 
the  Ith  set  of  creep  results 
are  to  be  printed 
( I MAX  =  NSPO ) 


I***:*::*:*:**::* 


X 


CARD  TYPE  NO.  14  -  CREEP  MATERIAL  PROPERTIES 

[Format(2(El5.4,F5.2) ,F5.2) ] 


Co  1 umns 

Variable 

1-15 

C0EF1 ( I ) 

16-20 

EXP  1  ( I  ) 

21-35 

C0EF2 ( I  ) 

36-40 

EXP2 ( 1  ) 

41-45 

EXP  3 (  1  ) 

-one  card  for  each  material 
-the  general  form  of  the  power 


Entry 

sf:  sf:  sfc  sfc 

Coefficient  in  the  creep  law 

Stress  exponent  in  the  creep 
1  aw 

Coefficient  in  the  creep  law 

Stress  exponent  in  the  creep 
1  aw 

Time  exponent  in  the  creep  law 

law  creep  relationship  is: 


e 

e 


ni  n  o 

[A,  a  +  A9a  ]  tm 
L  1  e  2  e  J 


where  e  =  effective  strain  rate 
e 

a  =  effective  strain 
e 

t  =  time 

A-j ,  A^  =  coefficients 
n-j,  n^  =  creep  exponents 
m  =  time  exponent 


. 


CARD  TYPE  NO.  15  -  CREEP  SOLUTION  CONVERGENCE  PARAMETERS 

[ Format ( 3  f 1 0 . 3 ) ] 


Columns  Variable 

*******  ******** 


Entry 


1-10  ETAO 


11-20  ETA  1 


21-30  OMEGA 


Maximum  ratio  of  effective 
elastic  strain  to  increment  of 
effective  creep  strain  for 
first  time  interval 

Maximum  ratio  of  effective 
elastic  strain  to  increment  of 
effective  creep  strain  for 
succeeding  time  intervals 

Maximum  fractional  change  in 
effective  stress  allowed  per 
creep  increment 


:  *  3fc  *  *  * 


-to  ensure  convergence,  ETAO,  ETA  1  and  OMEGA  should  lie  within 
the  following  ranges: 

10  <  ETAO  <  25 
1  <  ETA  1  <  10 
0.03  <  OMEGA  <0.10 

-default  values  have  been  set  equal  to: 

ETAO  =  25 
ETA  1  =  10 
OMEGA  =0.03 


-these  default  values  will  produce  the  minimum  time  increment 


X 


1  70 


CARD  TYPE  NC.  16  -  NODE  PRINT  OUT  CONTROL  (CREEP  SOLUTION  ONLY) 

[ Format (1216)] 


Columns  Variable 
*******  ******** 


Entry 


1-6  NNPR(I) 


7-12  NNPRU  +  1) 


Node  print  out  control 
indicator  for  node  number  I 
NNP R ( I ) =  0  print  out  and 

no  rotation 

NNPR(I)=-1  do  not  print 
NNPR ( I ) =  K  print  out  and 

rotate  an  angle 
of  K  degrees 
( k=an  integer ) 

Node  print  out  control 
indicator  for  node  number  1+1 


13-18  NNPR ( I +2 ) 


Node  print  out  control 
indicator  for  node  number  1+2 


67-72  NNPR ( 1+5 ) 


Node  print  out  control 
indicator  for  node  number  1+5 


:**************** 


-twelve  nodes  per  card 


. 
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CARD  TYPE  NO.  17  - 


ELEMENT  PRINT  OUT  CONTROL 
ONLY )  [ Format  (  1216) ] 


I******:*:*:*:***::*:*:*:** 

(CREEP  SOLUTION 


Columns  Variable 

*  *  *  sf:  *  *  * 


Entry 

:+:  :+c  3(e 


1-6  NEPR(I) 


7-12  NE  PR ( 1  +  1  ) 

13-18  NEPR ( 1+2 ) 


67-72  NEPR ( 1+5 ) 


Element  print 
indicator  for 
number  I 
NEPR (  I  )=  0 

NEPR ( I  )  =-  1 
NEPR ( I ) =  K 


Element  print 
indicator  for 
number  1+1 

Element  print 
indicator  for 
number  1+2 


Element  print 
indicator  for 
number  1+5 


out  control 
e 1 ement 

Print  out  and 
no  rotation 
do  not  print 
print  out  and 
rotate  an  angle 
of  K  degrees 
(k=an  integer) 

out  control 
e 1 ement 


out  control 
e 1 ement 


out  control 
e 1 ement 


-twelve  elements  per  card 


‘ 


APPENDIX  D 


PROGRAMME  LISTING  FOR  CREEP 
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APPENDIX  E 

THICK  WALL  CYLINDER  CLOSED  FORM  SOLUTION 


This  appendix  presents  a  comparison  of  the  finite  element  results 
predicted  by  the  programme  CREEP  and  the  closed  form  steady  state  creep  of  a 
thick  wall  cylinder  subjected  to  a  uniform  internal  pressure. 

The  finite  element  grid  is  shown  in  Figure  E— 1.  The  mesh  consists  of 
90  elements  and  61  nodes.  The  ratio  of  outer  to  inner  radius  is  egual  to 
2.0.  The  vertical  and  horizontal  boundaries  are  placed  on  ro I  I er  supports 
which  allow  translation  in  the  radial  direction  only. 

The  inner  radius  is  subjected  to  a  uniform  pressure  of  100  kPa.  The 
elastic  material  properties  are: 

E  =  10  x  106  kPa 
v  =  0.495 

The  flow  law  of  the  material  is; 

e  =  (1.0  x  10“^4)q  ^  ( e _ 1 ) 

_4  ^  5 

where  the  units  are  kPa  *  and  hr“1. 

The  steady  state  creep  velociy  for  the  inner  radius  is  shown  in  Figure 
^  2.  The  closed  form  solution  neglects  elastic  response  and  thus  maintains 
a  constant  value.  After  approximately  0.75  hr  of  creep  simulation,  the 
finite  element  solution  predicts  a  constant  velocity  2.2%  lower  than  the 
closed  form  solution. 

The  closed  form  and  finite  element  elastic  and  steady  state  effective 
stresses  are  compared  in  Figure  E-3.  As  shown  in  the  figure,  there  is  very 
good  agreement  for  both  elastic  and  steady  state  creep  stresses. 
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FIGURE  E-1 


THICK  WALL  CYLINDER  FINITE  ELEMENT  MESH 
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FIGURE  E-2  STEADY  -  STATE  CREEP  VELOCITY 


